Optimalizace algoritmu pro výpočet průniků simplexových výpočetních sítíOptimization of an algorithm for calculating the intersection of simplicial computing meshes

TECHNICKÁ UNIVERZITA V LIBERCI

Fakulta mechatroniky, informatiky a mezioborových studií

Studijní program: B2646 – Informační technologie Studijní obor: 1802T007 – Informační technologie

Optimalizace algoritmu pro výpočet průniků simplexových výpočetních sítí

Optimization of an algorithm for calculating the intersection of simplicial computing meshes

Bakalářská práce

Autor:

Viktor Friš

Vedoucí práce: Mgr. Jan Březina, Ph.D.

Vzhledem ke správnému číslování stránek nebude tato strana obsahem BP a bude nahrazena oficiálním zadáním.

Prohlášení

Byl(a) jsem seznámen(a) s tím, že na mou bakalářskou práci se plně vztahuje zákon č. 121/2000 Sb., o právu autorském, zejména § 60 – školní dílo.

Beru na vědomí, že Technická univerzita v Liberci (TUL) nezasahuje do mých autorských práv užitím mé bakalářské práce pro vnitřní potřebu TUL.

Užiji-li bakalářskou práci nebo poskytnu-li licenci k jejímu využití, jsem si vědom povinnosti informovat o této skutečnosti TUL; v tomto případě má TUL právo ode mne požadovat úhradu nákladů, které vynaložila na vytvoření díla, až do jejich skutečné výše.

Bakalářskou práci jsem vypracoval(a) samostatně s použitím uvedené literatury a na základě konzultací s vedoucím bakalářské práce a konzultantem.

Datum

Podpis

Poděkování

Rád bych poděkoval vedoucímu mé Bakalářské práce Mgr. Janu Březinovi, Ph.D. za poskytnutí odborných rad, ochotu, věcné připomínky, časté konzultace a vstřícný přístup během zpracování této práce.

Mé poděkování patří i mému dlouholetému kamarádovi a spolubydlícímu Jardovi Zlámalovi za podporu a motivaci práci dokončit.

Velké poděkování náleží celé mé rodině za podporu, trpělivost a povzbuzování po dobu mého studia.

Abstrakt

Bakalářská práce je zaměřena na optimalizaci algoritmu pro výpočet průniků simplexových výpočetních sítí oproti stávajícím algoritmům v programu Flow123d.

Sítě jsou tvořeny 3D elementy (čtyřstěny) a 1D elementy (úsečkami). Optimalizace algoritmu je dosaženo prodlužováním spočtených průniků a procházením sítě 3D elementů do šířky. Průniky jsou pro konkrétní dvojici 1D a 3D elementu vypočteny efektivním algoritmem využívající Plückerovy souřadnice. Další optimalizací je hledání kandidátů pro první průnik konkrétního jednoho elementu.

Pro efektivní výpočty průniků přímky s trojúhelníkem jsou zavedeny vztahy využívající Plückerovy souřadnice, které jsou následně dokázány. Průniky jsou reprezentovány lokálními souřadnicemi v souřadném systému 3D elementu nebo 1D elementu.

Pro práci se sítí se využívají datové struktury z programu Flow123d, zejména pro reprezentaci elementů a vrcholů, které se načítají z textového souboru.

Pro reprezentaci čtyřstěnů, úseček a dalších pomocných objektů se používají vlastní datové struktury.

Optimalizovaný algoritmus je detailně popsán a výsledky časové náročnosti oproti současnému algoritmu jsou dále prezentovány.

Implementace algoritmů a dalších funkcí jsou provedeny v jazyce C++.

Klíčová slova

Plückerovy souřadnice, průnik simplexů, nekompatibilní sítě, Flow123d

Abstract

The thesis is focused on optimization of the algorithm for computing the intersections of simplicial meshes over the algorithm currently used in the simulation software Flow123d. The meshes are made up of 3D elements (tetrahedrons) and 1D elements (lines). Optimization of the algorithm is achieved by extending already computed intersections and using the breadth-first search on the mesh of 3D elements. Intersections for specific pair of 1D and 3D elements are calculated by efficient algorithm using Plücker coordinates. Another optimization is the search for candidates for the first intersection of a particular element.

For the efficient calculation of the intersection of line and triangle are introduced relations using the Plücker coordinates. The relations are proven.

Intersections are represented by the local coordinates on the 3D element or by the local coordinates on the 1D element.

Data structures of the software Flow123d are used for work with meshes, especially for the representation of elements and vertices that are loaded from a text file. Own data structures are used to represent tetrahedra, lines, and other auxiliary objects.

The optimized algorithm is described in detail and the results of timings compared to the present algorithm are presented.

Algorithms and other functions are implemented in C++.

Keywords

Plücker coordinates, simplex intersections, non-compatible meshes, Flow123d

Obsah

Prohlášení...3

Poděkování...4

Abstrakt...5

Abstract...6

1 Úvod...9

2 Vlastnosti Plückerových souřadnic...11

2.1 Průnik přímky s trojúhelníkem, výpočet pomocí Plückerových souřadnic. .12 2.2 Invariance Plückerových souřadnic vzhledem ke změně počátečního bodu 13 2.3 Klasifikace vzájemné polohy přímek pomocí součinu Plückerových souřadnic...13

2.4 Důkaz výpočtu lokálních souřadnic pomocí Plückerových souřadnic...14

3 Datové struktury pro reprezentace sítí...17

3.1 Struktury programu Flow123d...17

3.2 Vlastní struktury...18

3.2.1 Simplex...18

3.3 Program GMSH...20

4 Lineární algoritmus pro hledání 1D-3D průniků...21

4.1 Prodlužovací bod...22

4.2 Algoritmus pro výpočet průniku 1D elementu se stěnou 3D elementu...23

4.3 Hlavní výpočetní třída...24

4.4 Hlavní cyklus průchodu přes všechny elementy...26

4.4.1 Inicializace obalových boxů pomocí BIHTree...26

4.4.2 Lineární inicializace obalových boxů...26

4.5 Inicializace fronty...27

4.6 Zpracování prodlužovacích bodů...28

4.7 Procházení 1D elementů do hloubky v jednom 3D elementu...31

5 Efektivita algoritmů...33

5.1 Diskuze o rychlosti algoritmů...35

6 Závěr...37

Seznam použité literatury...38

Obsah přiloženého CD...39

Seznam obrázků

Obr. 1: Možná relativní orientace dvou přímek p,q...11

Obr. 2: Normála k rovině A, B, B+V...13

Obr. 3: Simplex(S2)...15

Obr. 4: Simplex(S3)...19

Obr. 5: Příklad sítě v programu GMSH...20

Obr. 6: Vývojový diagram výpočtu prodlužovacího bodu...28

Obr. 7: Vývojový diagram zpracování prodlužovacího bodu...30

Obr. 8: Vývojový diagram zpracování 1D elementů uvnitř jednoho 3D elementu. .32 Obr. 9: Graf celkové časové náročnosti algoritmů pro různé sítě...34

Seznam tabulek

Tabulka 1: Číslování stěn a hran v simplexu S3...19

Tabulka 2: Celková časová náročnost algoritmů pro různé sítě...34

Tabulka 3: Časová náročnost inicializační části algoritmů...35

Tabulka 4: Časová náročnost výpočetní části algoritmů...35

1 Úvod

Puklinového proudění metodou konečných prvků je vytvářeno pro sítě různých dimenzí. Takové proudění je v tomto případě možné popsat pomocí komponent 1D elementů (úseček) v prostoru tvořeném 3D elementy (čtyřstěny).

Protože takové komponenty nemusí být pouze kompatibilní (takové, které proudí po hranách 3D elementů), ale mohou být i nekompatibilní (proudění prochází skrze 3D elementy), je potřeba průniky nekompatibilních komponent s 3D elementy nalézt, vypočítat a vhodně reprezentovat.

V této bakalářské práci se používá program Flow123d, jehož manuál lze nalézt na [1]. V programu je již implementován algoritmus pro hledání takových průniků. Algoritmus si vytváří binární strom sousedících elementů tzv. BIHTree, vysvětlen v článku [2], pro hledání kandidátů na průnik pro jeden konkrétní element (v našem případě 1D element). Časová složitost vytváření BIHTree je v ideálním případě O(N log N), kde N je celkový počet elementů v síti. Jednotlivé průniky pro konkrétní dvojici 1D elementu a 3D elementu se vypočítají pomocí soustavy rovnic přímky a trojúhelníku pro každou stěnu čtyřstěnu.

Cílem bakalářské práce je optimalizace algoritmu pro výpočet průniků simplexových výpočetních sítí oproti algoritmům implementovaných v programu Flow123d, z kterého se využívají nástroje pro práci s výpočetními sítěmi jako je jejich načítání a jejich reprezentace pomocí dalších datových struktur. Optimalizace algoritmu spočívá v prodlužování vypočtených průniků do sousedních elementů a procházení těchto elementů do šířky. Samotný výpočet průniků, pro konkrétní dvojici 1D a 3D elementu, je realizován efektivním výpočtem průniku přímky a trojúhelníku za použití Plückerových souřadnic pro každou stěnu čtyřstěnu.

Efektivní výpočet je inspirován z článku [3]. Druhou optimalizací je vytvoření jiného způsobu pro hledání kandidátů na průnik pro jeden element, než je použití BIHTree, jehož časová složitost by byla O(k N), kde N je celkový počet elementů v sítí a k je počet komponent 1D elementů v síti.

Plückerovy souřadnice jsou vysvětleny v kapitole 2, kde jsou odvozeny a dokázány jednotlivé vztahy pro efektivní výpočet průniků přímky s trojúhelníkem.

Dále jsou popsány jednotlivé datové struktury programu Flow123d, které jsou využity v algoritmech. V kapitole 3 jsou popsány vlastní vytvořené struktury pro práci s elementy ze sítě, zejména s čtyřstěnem. Jednotlivé optimalizace a algoritmus samotný je detailně popsán v kapitole 4, kde se řeší všechny důležité metody pro výpočet, hlavně procházení elementy do šířky a prodlužování průniků, jsou zde i popsány další struktury potřebné pro prodloužení a reprezentaci průniků.

Současný i optimalizovaný algoritmus bude otestován na pěti různých výpočetních sítích. Časové výsledky testů pro jednotlivé sítě jsou popsány v kapitole 5.

Implementace všech algoritmů, datových struktur a samotný program Flow123d je vytvářen v jazyce C++ [4].

2 Vlastnosti Plückerových souřadnic

Plückerovy souřadnice jsou speciální reprezentací přímky v třírozměrném prostoru za použití jistého šesti-rozměrného vektoru. Máme-li přímku p určenou počátečním bodem A a směrovým vektorem U, plückerova souřadnice přímky p je vypočtena tímto vztahem:

π p = (Up , Vp) = (U , U × A). (1)

Pro naše případy budeme potřebovat výpočet skalárního součinu Plückerových souřadnic dvou přímek p a q:

π p ʘ π q = Up ∙ Vq + Uq ∙ Vp , (2) který se liší oproti skalárnímu součinu vektorů, reprezentující dvě přímky (Up ∙ Uq + Vq ∙ Vp), v prohození vektorů Uq a Vq.

Výsledek indikuje orientaci přímky q vůči směru přímky p, pokud je orientace přímky p směrována do „hodin“, pak je přímka q ve vztahu podle Obr. 1.

Obr. 1: Možná relativní orientace dvou přímek p,q

(a) π p ʘ π q > 0 => q jde v protisměru hodinových ručiček kolem p.

(b) π p ʘ π q < 0 => q jde ve směru hodinových ručiček kolem p.

(c) π p ʘ π q = 0 => q protíná p nebo je s ní paralelní.

2.1 Průnik přímky s trojúhelníkem, výpočet pomocí Plückerových souřadnic

Součin dvou Plückerových souřadnic je důležitý při výpočtu průniku přímky a trojúhelníku. Pokud uvažujeme přímku p a trojúhelník dán vrcholy (V0,V1,V2) s orientací hran e0 = (V1,V2), e1 = (V2,V0), e2 = (V0,V1). Přímka p protíná trojúhelník, pokud skalární součiny Plückerových souřadnic přímky s každou hranou trojúhelníku jsou všechny buď kladné nebo záporné. Výhodou Plückerových souřadnic je i výpočet barycentrických souřadnic průsečíku v trojúhelníku, který je dán vztahem:

u_i^k=w_i^k/

∑

i=0 3

w_i^k , kde w_i^k=π_pʘ π_ei . (3) Zavedeme si lokální souřadnice trojúhelníku α a β, které budou reprezentovat průsečík v souřadném systému trojúhelníku. Průsečík v globálních (kartézských) souřadnicích vyjádříme obecným vztahem:

^Pk=V₀αe₂−βe₁ . (4)

Přičemž podle Tvrzení 3, odvodíme vztah pro výpočet globálních souřadnic průsečíku za použití Plückerových souřadnic a dostaneme vztah:

Pk=u0k

V0u1k

V1u2k

V2 . (5)

Pomocí lokálních souřadnic trojúhelníku lze zjistit lokální souřadnici průsečíku přímky. Uvažujme přímku danou bodem A a směrovým vektorem U, průsečík X, trojúhelník dán vrcholy (V0,V1,V2) a jeho lokální souřadnice α a β, potom platí: AθU = X =V₀αv₀−βv₂ , kde θ je lokální souřadnice přímky.

Po úpravě vznikne vztah: ^θ=V0−Aαv₀−βv₂/U. (6)

2.2 Invariance Plückerových souřadnic vzhledem ke změně počátečního bodu

Tvrzení 1: Uvažujme přímku r určenou bodem A a směrovým vektorem U.

Dále uvažujme tutéž přímku, ale určenou bodem X a směrovým vektorem U.

Plückerovy souřadnice budou stejné.

Důkaz: πr = (U , U × A) = (U , U × X) => U × A = U × X ? Bod X leží na stejné přímce, tedy X = A + αU. Proto po dosazení dostáváme

U × X = U × (A + αU) = U × A + α(U × U) = U × A + 0α = U × A a tedy (U , U × A) = (U , U × X).

Dokázáno

2.3 Klasifikace vzájemné polohy přímek pomocí součinu Plückerových souřadnic

Uvažujme přímku r určenou bodem A a směrovým vektorem U, jejíž plückerova souřadnice je πr = (U , U × A). Dále uvažujme přímku s určenou bodem B a směrovým vektorem V s plückerovou souřadnicí πs = (V , V × B). Jejich součin podle vztahu (2) je určen následovně:

π r ʘ π s = U ∙ (V × B) + V ∙ (U × A).

Součin Plückerových souřadnic znázorňuje směr obíhání jedné přímky kolem druhé nebo jejich protnutí, rovnoběžnost nebo shodu. To se pokusíme dokázat pomocí normály N k rovině A, B, B+V (viz Obr. 2). Normála N je určena vektorovým součinem směrového vektoru V a směrového vektoru BA, N = (A – B)×V.

Obr. 2: Normála k rovině A, B, B+V

Tvrzení 2: Skalární součin normály N se směrovým vektorem U bude odpovídat skalárnímu součinu dvou Plückerových souřadnic.

Důkaz: N ∙ U = π r ʘ π s ?

LS = ((A - B) × V) ∙ U = (A × V) ∙ U – (B × V) ∙ U = U ∙ (V × B) + (U × A) ∙ V PS = U ∙ (V × B) + V ∙ (U × A)

Dokázáno

Pokud se budou přímky protínat, normála bude kolmá na vektor U, tudíž skalární součin bude roven nule. Když budou přímky totožně, normála bude nulová a součin opět bude nulový a zároveň, když budou přímky rovnoběžné, tak bude normála kolmá k vektoru U a opět bude součin nulový. Bude-li součin větší než nula, normála s vektorem U svírá ostrý úhel, tudíž přímka s bude přímku r obíhat v protisměru hodinových ručiček. Naopak, když bude normála svírat tupý úhel bude přímka s obíhat přímku r ve směru hodinových ručiček.

N ∙ U = 0 => s, r jsou koplanární.

N ∙ U > 0 => s jde v protisměru hodinových ručiček okolo r.

N ∙ U < 0 => s jde ve směru hodinových ručiček okolo r.

Výsledky s obíháním jsou stejné jako u součinu Plückerových souřadnic a i jejich rovnost byla dokázána.

2.4 Důkaz výpočtu lokálních souřadnic pomocí Plückerových souřadnic

Tvrzení 3: Uvažujme přímku r určenou bodem X a směrovým vektorem U.

Dále uvažujme trojúhelník (viz Obr. 3), dán vrcholy (V0, V1, V2). Vztah pro barycentrické souřadnice průsečíku vzhledem k trojúhelníku bude záviset na skalárním součinu Plückerových souřadnic přímky s každou hranou trojúhelníku, z kterých se vypočítají lokální souřadnice.

Důkaz: u_i=π_rʘ π_vi/

∑

i=0 3

π_rʘ π_vi ? Kde ui jsou barycentrické souřadnice.

Dokážeme, že daný vztah platí. Počáteční bod přímky můžeme posunout do bodu průsečíku podle Tvrzení 1. Tudíž počáteční bod X na přímce je zároveň průsečíkem.

Průsečík vyjádříme vztahem:

X = V0 + (αv0 – βv2), (7)

kde lokální souřadnice α náleží barycentrické souřadnici u1 ve směru vektoru v0

a lokální souřadnice β náleží barycentrické souřadnici u2 ve směru vektoru v2.

Pro skalární součin Plückerových souřadnic přímky s každou hranou trojúhelníku dostáváme vztah:

w_i=π_rʘ π_vi=U v_i×V_iv_iU × X  . (8) Po dosazení vztahu (7) do vztahu (8) dostáváme rovnici pro první součin přímky r a hrany v0:

w₀=U v₀× V₀v₀U ×[V₀αv₀−βv₂]=U v₀×V₀v₀U × V₀

v₀U ×[αv₀−βv₂]=U v₀×V₀−U  v₀×V₀−U v₀×[αv₀−βv₂]=−βU  v₂× v₀.

Pro druhý součin přímky r a hrany v1 dostaneme rovnici:

Obr. 3: Simplex(S²)

w₁=U v₁×V₁v₁U ×[V₀αv₀−βv₂]=U v₁×[V₀v₀]v₁U ×V₀

v₁U ×[αv₀−βv₂]=U v₁× v₀v₁U ×[αv₀−βv₂].

vektor v1 vyjádříme pomocí ostatních dvou vektorů v0 a v1: v1 = - v0 – v2.

w₁=−U v₂×v₀−αU v₀× v₂βU  v₂× v₀=−U v₂× v₀α β [U v₂× v₀].

Pro třetí součin přímky r a hrany v2 dostaneme rovnici:

w₂=U  v₂× V₂v₂U ×[V₀αv₀−βv₂]=U v₂×V₀v₂U ×V₀

v₂U ×[αv₀−βv₂]=αU  v₀× v₂=−αU  v₂× v₀.

Všechny tři součiny jsme si upravili, aby byly ve stejném tvaru a lišily se jen konstantou.

w₀=−βU v₂×v₀ w₁=α β−1[U  v₂× v₀] w₂=−αU v₂×v₀

Když máme takto vyjádřené součiny, vyjádříme z nich lokální souřadnice α a β.

Při součtu všech tří rovnic dostaneme výsledek: ^w0w₁w₂=U v₁× v₂−1.

Z něhož již lehce dokážeme, že ^α=w2/w₀w₁w₂=u₁

^β=w0/w₀w₁w₂=u₂ (9) ^u0=1−u₁−u₂ .

Dokázali jsme že ui jsou barycentrické souřadnice a α a β jsou lokální souřadnice vypočteny ze skalárních součinů Plückerových souřadnic.

V závislosti na lokálních souřadnicích, odvodíme vztah pro průsečík, aby závisel pouze na vrcholech trojúhelníku a barycentrických souřadnicích.

X = V0 + (αv0 – βv2) = (1 – α – β)V0 + αV1+ βV2 = u0V0 + u1V1 + u2V2 .

3 Datové struktury pro reprezentace sítí

V této části se popíšou datové struktury pro popis sítí. Nejprve se popíšou struktury z programu Flow123d, pak struktury pro vlastní algoritmus. Bude popsán i program GMSH, který se využívá pro tvorbu a zobrazování sítí.

3.1 Struktury programu Flow123d

Významné struktury z programu Flow123d jsou zejména objekty pro práci se sítí. Jako jsou třídy Element, Edge, Side, Node, BIHTree, BoundingBox a Mesh:

Node – vrchol, který uchovává informaci o globální pozici vrcholu, reprezentovanou vektorem kartézských souřadnic. Je vytvářen načtením vrcholů ze sítě.

Element – Obsahuje zejména dimenzi daného elementu, jeho index v síti, počet vrcholů, odkazy na vrcholy (Node), počet stěn a konkrétní indexy stěn.

Je vytvářen načtením elementu ze sítě.

Side – stěna elementu. Dočasný objekt reprezentující jednotlivé stěny elementu. Obsahuje zejména odkaz na Edge.

Edge – objekt stěn. Stěny určené stejnými vrcholy nemusí být stejné objekty.

Každá stěna patří svému elementu. Pro rozlišení, které stěny jsou stejné vrcholy, slouží tato třída. Odkazuje na jednotlivé stěny sousedních elementů a tak lze skrze ni, zjistit sousední element. Obsahuje indexy na stěny jednotlivých sousedních elementů.

BoundingBox – obalový box nad jakýmkoliv elementem nebo množinou elementů. Je to nejmenší kvádr obalující element, jehož hrany jsou rovnoběžné se souřadných systémem (kartézským souřadných systémem).

BIHTree – strom obalových boxů vytvářen ze sítě, založen na jistém půlení intervalů. Dokáže zjistit indexy okolních obalových boxů, které se protínají s konkrétní obalovým boxem. BIHTree je vysvětlen v článku [2].

Mesh – Reprezentace celé sítě. Obsahuje veškerá data včetně vrcholů, elementů, stěn a dalších. Třída, která načte celý textový soubor se sítí z programu GMSH. .

3.2 Vlastní struktury

Pro realizaci efektivních výpočtů se používají operace s Plückerovými souřadnicemi, proto se navrhly třídy pro práci s nimi. Mezi základní třídy patří:

SPoint – reprezentuje bod v třírozměrném prostoru.

Vector – geometrický vektor, vytvořen z dvou bodů typu SPoint. Obsahuje metody pro výpočet skalárního a vektorového součinu.

Plucker – stěžejní třída, které dokáže pro libovolné dva body, tvořící přímku, spočítat Plückerovy souřadnice a uložit je. Obsahuje metodu pro výpočet součinu Plückerových souřadnic. Plückerovy souřadnice jsou reprezentovány šesti-rozměrným vektorem souřadnic.

HyperPlane – reprezentace přímky. Vytvořena z dvou bodů. Obsahuje i vlastní Plückerovy souřadnice.

Simplex – Třída reprezentující obecnou specifikaci n-rozměrného trojúhelníku (detailně vysvětleno níže).

3.2.1 Simplex

Simplex je reprezentací trojúhelníku v určitém prostoru. Pro naše konkrétní účely uvažujeme Simplex Sⁿ , kde n je dimenze, n = 0, 1, 2, 3. Hranice n-rozměrného simplexu jsou určeny (n-1)-rozměrnými simplexy. Pro představu simplex S⁰ je bod, simplex S¹ úsečka, simplex S² trojúhelník a simplex S³ čtyřstěn.

Každý simplex obsahuje n+1 vrcholů.

Třídu Simplex, schopnou vytvořit n-rozměrný simplex v určitém podprostoru, budeme zejména využívat pro tvorbu čtyřstěnů S³. Pro takový čtyřstěn se vytvoří simplexy stěn, hran a vrcholů. Libovolný simplex lze vytvořit předáním

pole vrcholů nebo nastavením (n-1)-rozměrnými simplexy. Uvažujme body v prostoru A, B, C a D, kterými určíme konkrétní simplex S³. Ten si vytvoří (n-1)-rozměrné simplexy s pevně danými indexy a orientací hran (viz Tabulka 1).

Stěny simplexu Vrcholy stěny Orientace hran stěny

S0 A, B, C 0(AB), 1(BC), 2(CA)

S1 A, B, D 0(AB), 3(BD), 4(DA)

S2 A, C, D 2(CA), 5(CD), 4(DA)

S3 B, C, D 1(BC), 5(CD), 3(BD)

Tabulka 1: Číslování stěn a hran v simplexu S³

Orientace hran je důležitá pro správné spočtení součinu Plückerových souřadnic a následné zjištění průniku ve stěně. Jelikož celá síť je načtena a zpracována do objektů typu Element a v programu Flow123d je záměr používat

Obr. 4: Simplex(S³)

síť založenou na podobné reprezentaci elementů jako je Simplex, je nutné zajistit správnou kompatibilitu mezi třídou Element a Simplex. Simplex se vytvoří načtením vrcholů ze třídy Element, bohužel číslování stěn je v opačném pořadí.

Pro Simplex je to 0, 1, 2 a 3 (viz Obr. 4), kdežto pro Element je to 3, 2, 1, 0.

Správná konverze je zajištěna přepočítáním indexů. Uvažujme index i stěny simplexu a index j stěny elementu, potom se indexy přepočítají vztahem: i = 3 – j.

3.3 Program GMSH

Grafický editor využíván pro tvorbu sítí (viz Obr. 5). Do programu se vkládají body v prostoru, které se následně mohou propojovat v úsečky a další geometrické entity. Z úseček se vytváří plochy a z ploch objemy. Entity se mohou ukládat do různých vrstev a tak lze oddělit např. síť plnou 3D objektů od sítě 1D objektů. Každému vytvořenému bodu lze nastavit diskretizační parametr, který udává jemnost dělení při síťování obrazce. Tří rozměrné síťování vyplní všechny objemy čtyřstěny a plochy trojúhelníky atd. Výsledná síť se uloží do textového soubor typu *.msh. Soubor obsahuje údaje o všech vrcholech, jejich indexy a souřadnice v prostoru. Obsahuje i informace o všech elementech, také jejich indexy, ale i údaj ke které vrstvě patří, údaje o materiálu, kolik vrcholů obsahuje a indexy jednotlivých vrcholů. Každý vrchol a element je uložen na samostatném řádku v souboru.

Obr. 5: Příklad sítě v programu GMSH

4 Lineární algoritmus pro hledání 1D-3D průniků

V této kapitole se zabýváme algoritmem pro výpočet průniků nekompatibilních komponent v sítích. Uvažujeme velkou prostorovou výpočetní síť tvořenou 3D elementy (čtyřstěny) a v ní síť 1D elementů (úseček), která může být tvořena z nesouvislých komponent. Algoritmus je založen na procházení sítě 3D a 1D elementů do šířky. Pokud již byla nalezena dvojice 1D a 3D elementů, ke které byl spočten průnik, využije se prodlužování průniku. Prodlužování se realizuje pro dva typy konců průniků:

• Průnik končí na stěně 3D elementu – 1D element protne 3D element ve stěně a pokračuje dál za něj. Pro takovýto průnik se vytvoří prodlužovací bod (viz kapitola 4.1), který se uloží do fronty. Je nutné najít přes společnou stěnu sousední 3D element, do kterého průnik pokračuje. Pokud se jedná o okrajový 3D element v sítí, nebude pro něj existovat sousední element a průnik se již neprodlužuje (nevytvoří se prodlužovací bod).

• Průnik končí na konci 1D elementu (uvnitř 3D elementu) – 1D element končí uvnitř 3D elementu. Prodlužování průniku se realizuje, tak že jsou zjištěny všechny sousední 1D elementy, které vychází z koncového bodu 1D elementu. Pokud není nalezen žádný soused, průnik se již neprodlužuje. Nalezenému sousedovi se spočte průnik se stejným 3D elementem. Sousedů může být více nebo nemusí být žádný. Pro každého souseda je okamžitě spočten průnik, dokud nekončí ve stěně 3D elementu, aby se vytvořil prodlužovací bod.

Sousedi jsou okamžitě zpracováváni, aby došlo k úspoře času při práci se stejným 3D elementem.

Inicializace prvotního průniku pro konkrétní 1D komponentu je realizována nalezením 1D elementu z komponenty, který nemusí být striktně na začátku komponenty. O každém 1D elementu je zaznamenána informace, jestli již byl

zkoumám. Pro nalezený 1D element se hledá 3D element buď přes celou síť (viz kapitola 4.4.2), pro který bude existovat neprázdný průnik nebo pomocí BIHTree (viz kapitola 4.4.1).

4.1 Prodlužovací bod

Procházení sítě 3D elementů do šířky je realizováno díky pomocné datové struktuře nazvané prodlužovací bod. Poté, co se vytvoří průnik končící na stěně 3D elementu se vytvoří i prodlužovací bod, který patří k sousednímu 3D elementu.

Pokud sousední element neexistuje, prodlužovací bod se nevytvoří. Postupným prodlužováním průniků se vytváří prodlužovací body, které se ukládají do fronty.

Z fronty se pak načítají a slouží k výpočtu dalších průniků. Pokud leží celá 1D komponenta v 3D síti, bude se fronta procházet, dokud se nenaleznou průniky pro všechny elementy 1D komponenty. Do prodlužovacího bodu se ukládají informace o vzniklém průsečíku na stěně 3D elementu tak i o tom samém průsečíku vzniklém na 1D elementu. Ukládají se následující data:

• Indexy dvojice elementů – číselný index 1D a 3D elementu v síti, ke kterým prodlužovací bod patří. Stačí si ukládat pouze číselný index na oba příslušné elementy, z kterých se snadno určí konkrétní element.

• Číslo stěny – lokální číslo stěny 3D elementu. Pomocí něhož lze určit sousední 3D element. Navíc slouží i k optimalizaci, protože pro následný výpočet průniku nového 3D elementu se již nemusí procházet tato stěna.

• Orientace 1D elementu – pomocí výpočtu průniku 1D elementu se stěnou 3D elementu Plückerovými souřadnicemi (viz kapitola 4.2), lze určit orientaci 1D elementu, zda-li element vstupuje do 3D elementu nebo vystupuje. Pro naše účely si pevně určíme, že pokud vstupuje, uložíme si pravdivostní hodnotu true a pokud vystupuje, uložíme si pravdivostní hodnotu false.

• Lokální souřadnice průsečíku pro stěnu 3D elementu – souřadnice jsou uloženy do pole číselných hodnot v souřadném systému 3D elementu. Hodnoty jsou v intervalu <0;1>. Jelikož většina průsečíků vzniká na stěnách, stačí uchovávat pouze dvě hodnoty, pro zbylé přeinterpolované průsečíky se ukládají hodnoty tři.

• Lokální souřadnice průsečíku pro 1D element – pouze jediná číselná hodnota reprezentující pozici průsečíku na 1D elementu.

4.2 Algoritmus pro výpočet průniku 1D elementu se stěnou 3D elementu

Výpočet průniku pro konkrétní dvojici 1D a 3D elementu je realizováno převodem 1D elementu na objekt přímky a 3D elemenu na objekt čtyřstěnu S³. Pro tyto objekty jsou implementovány metody na výpočet průniku. Průnik se stěnou se spočte za využití součinu Plückerových souřadnic přímky s každou hranou stěny (viz kapitola 2.1). Pokud budou všechny tři součiny se stejnými znaménky, přímka protne stěnu. Když přímka skrze stěnu vstupuje do čtyřstěnu, budou součiny kladné.

Vystupuje, pokud součiny budou záporné. Jednotlivé hrany v jedné stěně se orientují všechny ve stejném směru. Vypočtené součiny se ukládají podle orientace hran v čtyřstěnu (viz Tabulka 1), proto je nutné orientace přepočítat. Přepočet orientace se realizuje přenásobením součinu mínus jedničkou. Výsledné tři součiny stěny, se načítají z pole součinů hran pro konkrétní čtyřstěn. Při načítání je opět nutné kontrolovat správnou orientaci hran, nejen aby byly všechny hrany ve stejném směru, ale aby hodnoty kladné odpovídaly pro vstup úsečky do čtyřstěnu a záporné pro výstup. Pro každou stěnu, ke které se hledá průnik, se počítají součiny Plückerových souřadnic, aby se nemuselo zbytečně počítat 12 součinů pro celý čtyřstěn, ukládají se součiny do pole o 6 hodnotách (viz kapitola 4.3), kde každá patří k hraně čtyřstěnu a případně se součin přenásobí, aby odpovídal orientaci v čtyřstěnu. Při následném načítání se kontroluje, zda-li součin už neexistuje a jestli ano, jestli je ve správné orientaci pro konkrétní stěnu. Ze součinů se vypočítají lokální souřadnice pro čtyřstěn i pro úsečku (viz vztah (6) a (9)).

4.3 Hlavní výpočetní třída

Veškeré výpočetní operace probíhají ve třídě InspectElements, která obsahuje metody pro výpočet průniků, prodloužení průniků zpracováním fronty prodlužovacích bodů, převod elementu na vlastní struktury atd. Třída obsahuje následující data:

• Pole pravdivostních hodnot prozkoumaných elementů – když se vytvoří síť s určitým počtem elementů, vytvoří se i stejně velké pole pravdivostních hodnot, které udává informaci o tom, zda-li byl již element zkoumán. Pro naše účely postačí jen tato informace o 1D elementech. Dříve než se vyšetřuje 1D element, jestli neprotíná 3D element nebo jestli se neprodlužuje průnik a hledá se nový 1D element, ověří se jeho index v tomto poli. Při cyklických topologiích 1D elementů v komponentách bude zaručeno, že se algoritmus nezacyklí. Pokud se 1D komponenta různě větví, je díky této informaci zaručeno, že se nebudou již procházet 1D elementy, které se prošly. Při hledání prvotního elementu dané komponenty se vyloučí elementy, pro které již byl průnik spočten.

• Pole průniků – výsledné dynamické pole, které obsahuje všechny spočtené průniky. Každý průnik je reprezentován třídou IntersectionLocal. Třída IntersectionLocal obsahuje indexy dvojic 1D a 3D elementu, ke kterým patří. Obsahuje i ukazatele na dva průsečíky z kterých je vytvořen průnik, tyto ukazatele jsou typu IntersectionPoint, pomocná třída, která přebírá lokální souřadnice pro 3D i pro 1D element. Při prodlužování průniku přes stěnu 3D elementu se všechny lokální souřadnice předávají z prodlužovacího bodu.

• Frontu – nezbytná fronta pro průchod 3D elementů do šířky. Ukládá prodlužovací body. Fronta je typu FIFO, prodlužovací bod, který dřív bude uložen, bude dříve zpracován.

• Odkaz na síť – vstupní síť je předána jako odkaz. Pomocí indexů z ní lze získat dané elementy, stěny, hrany atd. Obsahuje i některá makra na procházení všech jejich elementů.

• Pomocné data Simplex a HyperPlane – jelikož se v síti pracuje s objekty typu Element a ve výpočetních algoritmech se pracuje s objekty Simplex a Hyperplane, je nutný převod objektu Element na typ Simplex nebo HyperPlane. Simplex se vytvoří podle vrcholů 3D elementu vždy, když se hledá průnik pro novou dvojici 1D a 3D elementu. HyperPlane, přímka, která se vytváří z vrcholů 1D elementu. U přímky je nutností sledovat její správnou orientaci, proto při vytváření objektu se kontroluje informace o orientaci z prodlužovacího bodu. Pokud je orientace opačná, přímka se vytvoří s prohozenými body. Pokud při vytváření přímky není prodlužovací bod znám, orientace zůstane taková, jakou ji má 1D element.

• Pole součinů – pole o šesti hodnotách, kde každá odpovídá součinu Plückerových souřadnic 1D elementu s jednou hranou 3D elementu.

Součiny odpovídají hranám v čtyřstěnu v dané orientaci (viz Obr. 4).

Pole je dočasné jen pro práci s konkrétní dvojicí 1D a 3D elementu.

Při změně jakéhokoliv elementu jsou hodnoty pole nastaveny na hodnotu Null.

Třída InspectElements obsahuje navíc metody pro naplnění pole součinů, metody na interpolaci lokálních souřadnic ve 3D elementu a základní metody na získání polí a důležitou metodu pro výpis průniků do souboru a jejich následný převod do globálních souřadnic.

Jednotlivé metody k nalezení průniků a k prodlužování průniků jsou detailně vysvětleny níže.

4.4 Hlavní cyklus průchodu přes všechny elementy

Načtením souboru se sítí se vytvoří veškeré elementy, které je potřeba procházet a hledat první prodlužovací bod a pak je zpracovávat. Nad každým elementem se vytváří obalový box, díky kterému lze snadno zjistit kandidáta na průnik. Inicializace obalových boxů je realizována dvěma způsoby, které jsou pak časově testovány a jejich efektivita je promítnuta ve výsledcích (viz kapitola 5).

4.4.1 Inicializace obalových boxů pomocí BIHTree

Vytvoří se strom obalových boxů nad elementy pomocí datové struktury BIHTree z programu Flow123d, která umožní zjištění všech obalových boxů, které protínají konkrétní jeden obalový box. Po inicializaci obalových boxů se procházejí všechny elementy sítě a hledá se 1D element, který ještě nebyl prozkoumán z pole pravdivostních hodnot z hlavní třídy. Z nalezeného 1D elementu se vytvoří jeho obalový box, pro který se naleznou všechny indexy kandidátů obalových boxů na průnik ze stromu. Pole těchto indexů se prochází. Index obalového boxu je i indexem jeho elementu, proto se danému elementu zjistí dimenze a je-li element 3D, vypočte se prodlužovací bod (viz kapitola 4.5). Pokud pro daný 3D element byl spočten prodlužovací bod, prochází se fronta prodlužovacích bodů, dokud není prázdná. Z fronty se vybere prodlužovací bod, který do ní první vstoupil a zpracuje se (viz kapitola 4.6), po zpracování se prodlužovací bod odstraní z fronty.

4.4.2 Lineární inicializace obalových boxů

Vytvoří se pomocné pole obalových boxů nad všemi elementy. Během vytváření obalových boxů se vytváří i obalový box nad celou sítí 3D elementů.

Po inicializaci obalových boxů se procházejí všechny elementy sítě a hledá se 1D element, který ještě nebyl prozkoumán z pole pravdivostních hodnot z hlavní třídy.

Spočte se průnik obalových boxů 1D elementu s obalovým boxem celé 3D sítě, zda-li vůbec 1D element leží v sítí. Jelikož neexistuje žádný strom, který by zjistil kandidáty obalových boxů na průnik, musí se procházet znovu celá síť elementů.

Pro každý 3D element je zjištěn průnik jeho obalového boxu s obalovým boxem 1D elementu. Pokud je průnik neprázdný, vypočte se prodlužovací bod (viz kapitola 4.5). Pokud pro daný 3D element byl spočten prodlužovací bod, prochází se fronta prodlužovacích bodů, dokud není prázdná. Z fronty se vybere prodlužovací bod, který do ní první vstoupil a zpracuje se (viz kapitola 4.6), po zpracování se prodlužovací bod odstraní z fronty.

4.5 Inicializace fronty

Stěžejní část algoritmu je prodlužování průniků ze stěny 3D elementu a procházení sousedním 3D elementem. Pro realizace řešení je nutné vypočítat prodlužovací bod, který určuje konkrétní dvojici 1D a 3D elementu. Přes první spočtený prodlužovací bod lze zaručit výpočet průniků pro celou souvislou komponentu 1D elementů. Inicializace fronty je znázorněna na Obr. 6.

Metoda nejdříve převezme 1D a 3D element, pro který se hledá průnik a prodlužovací bod. Tyto elementy se musí převést na jiné objekty, pro které jsou implementovány metody pro výpočet průniků. 3D element se převede na objekt typu Simplex (dále jen čtyřstěn). 1D element se převede na objekt typu HyperPlane (dále jen úsečka). Jelikož nebyl spočten zatím žádný prodlužovací bod, nelze zjistit, jestli je úsečka v orientaci, která je potřeba na to, aby do čtyřstěnu vstupovala.

Proto, když se úsečka vytváří, zůstane ji orientace původní (ve směru prvního a druhého vrcholu). Vytvoří se pomocná dočasná číselná data, která budou naplněna informacemi o lokálních souřadnicích a o správné orientaci. Prochází se každá stěna čtyřstěnu, kde pro každou se vypočítává průnik (viz kapitola 4.2). Pokud je průnik neprázdný, naplní se pomocná data. Vzniklý průnik nemusí ležet na úsečce, takový průnik se dále nezpracovává, proto se porovná vypočtená lokální souřadnice úsečky, zda-li je v intervalu <0;1>. Pokud není, metoda vrátí hodnotu false, protože se nepodařilo vypočítat průnik a vytvořit prodlužovací bod. Metoda vrátí hodnotu false i pro případ, že velmi malý 1D element leží celý uvnitř 3D elementu. Pokud ovšem průnik spadá do intervalu, vytvoří se prodlužovací bod, který je naplněn hodnotami

z pomocných dat. Je mu i přidělena spočtená orientace, jestli úsečka byla správně orientovaná, když vstupovala do čtyřstěnu. Prodlužovací bod se uloží do fronty pro budoucí zpracování.

Jelikož se jedná o první prodlužovací bod na stěně 3D elementu, bude existovat prodlužovací bod i na stěně sousedního 3D elementu (pokud existuje sousední 3D element). Index stěny čtyřstěnu, ke kterému byl prodlužovací bod nalezen, je ale jiný než index stěny 3D elementu (viz kapitola 3.2.1). Po správném převodu stěn se nalezne sousední 3D element. Pokud existuje, vytvoří se druhý prodlužovací bod, naplní se stejnými hodnotami, pouze se změní orientace úsečky vzhledem k prvnímu prodlužovacímu bodu. Pokud úsečka vstupuje do čtyřstěnu, musí ze sousedního čtyřstěnu vystupovat a obráceně. Do druhého prodlužovacího bodu se musí změnit index stěny, protože přepočítaný index odpovídá stěně 3D elementu, nikoliv stěně čtyřstěnu. Prodlužovací bod se uloží do fronty. Metoda vrátí hodnotu true, pokud byl uložen jeden nebo dva prodlužovací body.

4.6 Zpracování prodlužovacích bodů

Prodlužovací bod obsahuje veškeré potřebné informace ke spočtení průniku pro danou dvojici 1D a 3D elementu. Pomocí něj se nemusí pro jeden 1D element procházet celá síť 3D elementů. Jeho zpracováním může vzniknout další Obr. 6: Vývojový diagram výpočtu prodlužovacího bodu

prodlužovací bod sousedního 3D elementu a tak je zaručeno, že se spočtou průniky pro celou souvislou komponentu 1D elementů. Zpracování prodlužovacích bodů je znázorněna na Obr. 7.

Metoda přebírá daný prodlužovací bod a zpracovává jej. Načte z něj indexy na konkrétní dvojici 1D a 3D elementů a pomocí nich načte ze sítě objekty těchto elementů. 3D element se převede na objekt čtyřstěnu a 1D element na objekt úsečky. Z prodlužovacího bodu je načtena informace o tom, jestli byla úsečka ve správné orientaci, podle toho se případně úsečka přeorientuje prohozením bodů.

Vytvoří se pomocná číselná data, která se naplní hodnotami z výpočtu průniku úsečky se stěnou čtyřstěnu. V prodlužovacím bodě jsou uloženy jedny lokální souřadnice z průniku úsečky se stěnou čtyřstěnu, ale pro vytvoření průniku úsečky se čtyřstěnem je potřeba znát dva takovéto průniky. Proto se procházejí všechny stěny čtyřstěnu, kromě té, ke které byl již jeden průnik spočten. Informaci o stěně, která už byla použitá, se zjistí z prodlužovacího bodu. Pro každou stěnu se vypočítá průnik, je-li neprázdný, naplní se pomocná data informacemi o průniku a poli pravdivostních hodnot v hlavní třídě (viz kapitola 4.3) se indexu 1D elementu nastaví hodnota true, že byl již prozkoumán. Vzniklý průnik nemusí ležet na úsečce, mohou nastat dva případy:

• Průnik leží na úsečce – lokální souřadnice průniku na úsečce je v intervalu <0;1>. Vytvoří se objekt průniku 1D-3D elementu, do kterého se uloží indexy obou elementů, lokální souřadnice průniku z prodlužovacího bodu a lokální souřadnice druhého průniku, právě spočteného. Pokud byla úsečka přeorientována, lokální souřadnice průniku na úsečce se musí přepočítat. Výsledný objekt průniku 1D-3D se uloží do pole všech nalezených průniků. Pokud existuje přes stěnu sousední 3D element, vytvoří se pro něj prodlužovací bod. Stěny 3D elementu se musí opět přepočítat na stěnu čtyřstěnu a obráceně.

Prodlužovací bod je naplněn hodnotami z vypočítaného průniku a uloží se do fronty. Informaci o orientaci úsečky se do nového

prodlužovacího bodu uloží ze zpracovávaného bodu.

• Průnik leží mimo koncové body úsečky – Jelikož úsečka je vždy orientována, že vchází do čtyřstěnu, bude lokální souřadnice úsečky vždy větší než jedna, protože se úsečka prodlouží ve svém směru a průnik tedy vznikne „před ní“. Taková to úsečka končí uvnitř 3D elementu, proto lokální souřadnice 3D elementu z průniku se musí přeinterpolovat, aby odpovídaly souřadnicím uvnitř čtyřstěnu, nikoliv na stěně. Podle orientace z prodlužovacího bodu bude lokální souřadnice na úsečce rovna nule nebo jedné. Vytvoří se objekt průniku 1D-3D elementu, do kterého se uloží indexy obou elementů, lokální souřadnice průniku z prodlužovacího bodu, lokální přeinterpolované souřadnice 3D elementu a lokální souřadnice 1D elementu (nula nebo jednička). Výsledný objekt průniku 1D-3D se uloží do pole všech nalezených průniků.

Pokud úsečka končí uvnitř čtyřstěnu, naleznou se všechny sousední 1D elementy přes koncový/počáteční bod. Indexy všech těchto 1D elementů jsou ověřeny v poli pravdivostních hodnot v hlavní třídě, zda-li již element nebyl zkoumán. Pokud nebyl, je každý zpracováván metodou níže.

Obr. 7: Vývojový diagram zpracování prodlužovacího bodu

4.7 Procházení 1D elementů do hloubky v jednom 3D elementu

Když uvnitř 3D elementu končí 1D element, může se komponenta větvit.

Tyto větve 1D elementů jsou okamžitě zpracovány a jsou vypočteny průniky pro stejný 3D element, dokud se nenalezne průnik na stěně 3D elementu a vytvoří se prodlužovací bod. Metoda je znázorněna na Obr. 8.

Metoda nepracuje s žádným prodlužovacím bodem a proto se ji musí spoustu parametrů předat. Předává se ji index na nový a předchozí 1D element, index na 3D element, přeinterpolované lokální souřadnice 3D elementu a orientace předchozího 1D elementu. Z indexů 1D elementů se vytvoří elementy. Musí se zajistit správná orientace nového 1D elementu, proto se porovnává předchozí orientace 1D elementu:

• Pokud byl 1D element ve směru vstupu do 3D elementu – koncový bod 1D elementu je uvnitř 3D elementu, pokud tento bod bude stejný jako počáteční bod nového 1D elementu, orientace bude v pořádku, jinak se 1D element musí přeorientovat.

• Pokud byl 1D element ve směru výstupu z 3D elementu – počáteční bod 1D elementu je uvnitř 3D elementu, pokud tento bod bude stejný jako počáteční bod nového 1D elementu, orientace bude v pořádku, jinak se 1D element musí přeorientovat.

Vytvoří se správně orientovaný objekt úsečky, objekt čtyřstěnu zůstane z předchozí metody nezměněn. Procházejí se všechny stěny čtyřstěnu, dokud se nenalezne neprázdný průnik. Pro takový se ověří orientace. Pokud jsou všechny úsečky správně orientovány, jediná validní orientace je taková, kdy úsečka míří ven z čtyřstěnu a tedy lokální souřadnice průniku na úsečce bude vždy větší než nula.

Pro takový průnik je v metodě uloženo číslo stěny a přes index 1D elementu se v poli pravdivostních hodnot v hlavní třídě zaznamená, že byl již element prozkoumán. Lokální souřadnice průniku na úsečce je vždy větší než nula, ale může být větší nebo menší než jedna:

• Souřadnice je menší než jedna – průnik leží na úsečce a byl vytvořen na stěně 3D elementu. Vytvoří se objekt průniku 1D-3D elementu, který se naplní indexy elementů a uloží se do něj předaná dvojice přeinterpolovaných lokálních souřadnic 3D elementu s lokální souřadnicí počátečního bodu na 1D elementu. Druhá dvojice lokálních souřadnic se předá z vypočítaného průniku na stěně a objekt průniku 1D-3D elementu se uloží do pole všech průniků. Pokud existuje přes stěnu sousední 3D element, vytvoří se pro něj prodlužovací bod, který se uloží do fronty.

• Souřadnice je větší než jedna – úsečka začíná i končí uvnitř 3D elementu. Lokální souřadnice na úsečce je nula nebo jedna. Lokální souřadnice na stěně 3D elementu se musejí přeinterpolovat. Vytvoří se objekt průniku 1D-3D elementu, který se naplní indexy elementů a uloží se do něj obě přeinterpolované lokální souřadnice v 3D elementu a obě lokální souřadnice na 1D elementu (hodnoty nula a jedna). Jelikož úsečka končí uvnitř čtyřstěnu, naleznou se všechny sousední 1D elementy přes koncový či počáteční bod. Indexy všech těchto 1D elementů jsou ověřeny v poli pravdivostních hodnot v hlavní třídě, zda-li již element nebyl zkoumán. Pokud nebyl, je každý zpracováván rekurzivně stejnou metodou.

Obr. 8: Vývojový diagram zpracování 1D elementů uvnitř jednoho 3D elementu

5 Efektivita algoritmů

Vytvořené sítě byly testovány pro 3 algoritmy. První algoritmus využívá strom BIHTree pro hledání kandidátů na průnik pro konkrétní jeden element.

Průniky vypočítává již implementovanou funkcí getIntersection() z programu Flow123d. Tato funkce nevyužívá žádné prodloužení průniků, pro každý 1D element najde kandidáty 3D elementů a průnik počítá pomocí soustav rovnic rovin a přímky. Druhý je lineární algoritmus pro hledání 1D-3D průniků (viz kapitola 4) za využití stromu BIHTree (viz kapitola 4.4.1). Podle průběžných časových testů (viz tabulka 2) nebyl algoritmus o moc efektivnější, proto podle časů vytváření stromu BIHTree (viz tabulka 3) byl zaveden třetí algoritmus, který místo použití stromu vytváří obalové boxy nad všemi elementy v síti a pro jednotlivé 1D elementy jedné komponenty hledá kandidáty na průnik přes celou síť (viz kapitola 4.4.2).

Byly testované sítě s různou topologií a různým počtem komponent 1D elementů. Každá síť byla testována několikrát pro každý algoritmus a průměrné výsledky celkových časů jsou promítnuty v tabulce (viz tabulka 2) a zobrazeny v grafu (viz Obr. 9).

Byly vytvořeny a použity následující sítě:

1. Síť – obsahuje dvě rozdílné komponenty 1D elementů. První je v trojúhelníkové topologii, které je z části uvnitř 3D sítě. Druhá je hvězdicové topologie, rozdělená na 4 větve, které zasahují i mimo 3D síť. 1D elementů je 10 % z celkového počtu elementů.

2. Síť – obsahuje dvě rozdílné komponenty 1D elementů. První je v trojúhelníkové topologii, která je celá uvnitř 3D sítě. Druhá je nerozvětvená cesta skrze 3D síť. 1D elementů je 8 % z celkového počtu elementů.

3. Síť – obsahuje jednu velkou komponentu 1D elementů v hvězdicové topologii, rozdělenou na 5 větví, která je uvnitř 3D sítě a jednotlivé větve

zasahují i mimo 3D síť. 1D elementů je 30 % z celkového počtu elementů.

4. Síť – obsahuje přes deset různých komponent 1D elementů. 1D elementů je 5 % z celkového počtu elementů. Komponenty zasahují dovnitř 3D sítě i ven.

5. Síť – obsahuje přes třicet různých komponent 1D elementů. 1D elementů je 5 % z celkového počtu elementů. Komponenty jsou některé celé mimo 3D síť, některé uvnitř a některé uvnitř i mimo síť.

Číslo sítě Počet elementů 1. algoritmus 2. algoritmus 3. algoritmus

1 164188 4,35 s 3,18 s 0,86 s

2 545806 17,05 s 11,18 s 2,51 s

3 701320 33,05 s 14,46 s 2,78 s

4 528438 22,52 s 19,87 s 9,16 s

5 527360 21,3 s 18,18 s 4,69 s

Tabulka 2: Celková časová náročnost algoritmů pro různé sítě

Efektivnost různých části algoritmů je znázorněna v tabulce 3 a 4. V první tabulce je znázorněna časová náročnost inicializační části, kde se vytváří buď strom BIHTree nebo se vytváří obalové boxy všech elementů (lineární inicializace).

V druhé tabulce je zobrazena časová náročnost výpočetní části algoritmu, která obsahuje procházení všemi elementy a hledání kandidátů na průnik pro 1D element i výpočet samotných průniků.

Obr. 9: Graf celkové časové náročnosti algoritmů pro různé sítě

1. Síť 2. Síť 3. Síť 4. Síť 5. Síť

0 5000 10000 15000 20000 25000 30000 35000

1. Algoritmus 2. Algoritmus 3. Algoritmus

ms

Číslo sítě Vytvoření BIHTree Lineární inicializace obalových boxů

1 2,7 s 0,3 s

2 10,64 s 1 s

3 11,35 s 1,1 s

4 10,11 s 0,99 s

5 9,75 s 0,98 s

Tabulka 3: Časová náročnost inicializační části algoritmů

Číslo sítě 1. algoritmus 2. algoritmus 3. algoritmus

1 1,65 s 0,48 s 0,56 s

2 6,41 s 0,54 s 1,51 s

3 21,7 s 3,11 s 1,68 s

4 12,41 s 8,76 s 8,17 s

5 11,55 s 8,43 s 3,71 s

Tabulka 4: Časová náročnost výpočetní části algoritmů

5.1 Diskuze o rychlosti algoritmů

Silnou stránkou 2. a 3. algoritmu je lineární výpočet průniků pro jednu komponentu 1D elementů. Čím méně je v sítí komponent a čím více elementů jednotlivé komponenty obsahují, tím budou oba algoritmy efektivnější. Velký časový nárok má zejména BIHTree strom, který se vytváří pro 1. a 2. algoritmus.

Jeho vytváření přináší výhodu jen, když síť obsahuje velké množství komponent.

Algoritmy byly testovány pro sítě, kde množství 3D elementů převažuje nad množstvím 1D elementů a počty komponent jsou v řádech desítek, protože předpokládáme, že tak to bude i pro reálné sítě (např. pukliny v horninách, toky horninou atd.).

Časy 3. algoritmu z tabulky 4 jsou včetně hledání kandidátů a výpočtů průniků dohromady, proto někdy převyšují časy 2. algoritmu. Lze si povšimnout, že pro 5. síť, kde je více komponent je výpočet rychlejší oproti 4. síti, kde je méně komponent, při podobných počtech elementů. Je to způsobeni tím, že i když je v 5. síti více komponent, některé sou celé mimo síť a pro jednotlivé 1D elementy z takové komponenty se nemusí hledat kandidáti.

Jak je vidět také z tabulky 4, tak vyhledávání kandidátů pro větší sítě už není tak výhodné jako je u 2. algoritmu. Hledání kandidátů přes všechny elementy je v rámci 3. algoritmu pro větší sítě výhodnější. Časová náročnost nalezení kandidátu je u 3. algoritmů dána počtem komponent, protože se hledají kandidáti přes všechny elementy. Dobře to je vidět u 4. sítě, kde je mnoho komponent uvnitř. Na rozdíl u 5. sítě je spousta komponent i mimo celou síť, tudíž pro takové komponenty se kandidáti nehledají a časová náročnost je téměř nulová.

Z tabulky 3 je patrné, že vytváření stromu BIHTree, oproti zbytků výpočtů v 2. algoritmu, je příliš časově náročné a pro podobné sítě, jako byly testované, se nevyplatí BIHTree vytvářet.

6 Závěr

Cílem práce bylo vytvoření optimalizovaného algoritmu pro výpočet průniků simplexových výpočetních sítí. Algoritmus byl optimalizován oproti algoritmu z programu Flow123D těmito způsoby:

1. Místo soustav rovnic, používaných v současném algoritmu pro výpočet průniku, byly použity Plückerovy souřadnice, které jsou podle článku [3]

efektivnější. Jejich vztahy pro výpočet průniků byly dokázány a byly implementovány do optimalizovaného algortimu.

2. Místo vytváření binárního stromu BIHTree a jeho použití pro hledání kandidátů na průnik pro jeden 1D element, kde pro každého kandidáta byl spočten průnik (časová složitost O(N log N), kde N je počet elementů v síti), byl použit algoritmus, který prochází síť sousedících 3D elementů do šířky a prodlužuje spočtené průniky (časová složitost O(k N), kde k je počet komponent 1D elementů).

Vytvořený algoritmus umožňuje procházet síť sousedících 3D elementů do šířky a v rámci jednoho 3D elementů procházet síť 1D elementů do hloubky, pro které následně spočte průniky. Algoritmus dokáže i vypočítat lokální souřadnice průniků na elementech a tím kompletně reprezentovat průnik oproti současnému algoritmu, který lokální souřadnice vypočítat nedokáže.

Pro algoritmus bylo vytvořeno 5 různých sítí, které obsahovaly řádově sta tisíce elementů, pro všechny byl optimalizovaný algoritmus efektivnější, v jistých případech byl rychlejší až o 91 % oproti stávajícímu algoritmu.

V optimalizování algoritmu lze pokračovat. Mohly by se předávat vypočtené Plückerovy souřadnice mezi sousedními 3D elementy nebo by se mohly Plückerovy souřadnice předpočítat pro každou hranu elementu pro celou síť atd.

Výpočetní části a datové struktury algoritmu by mohly být použity na vývoj lineárního algoritmu pro hledání 2D-3D průniků.

Seznam použité literatury

[1] BŘEZINA, Jan. Flow123D: Documentation of file formats and brief user manual [online]. Verze 1.6.5. [cit. 2011-10-25].

URL: <https://dev.nti.tul.cz/~brezina/flow_doc/flow123d_manual.pdf>.

[2] Wikipedia contributors. Bounding interval hierarchy [online].

[cit. 2013-05-14].

URL: <http://en.wikipedia.org/wiki/Bounding_interval_hierarchy>.

[3] Nikos Platis, Theoharis Theoharis. Fast Ray-Tetrahedron Intersection Using Plücker Coordinates [online]. [cit. 2003].

URL: <http://users.uop.gr/~nplatis/files/PlatisTheoharisRayTetra.pdf>.

[4] Learn C++. The C++ Tutorial [online]. [cit. 2007-05-25].

URL: <http://www.learncpp.com>.

Obsah přiloženého CD

Na přiloženém disku se nachází:

• Použité sítě

• Složka obsahující zdrojové kódy optimalizovaného algoritmu a vlastních datových struktur

• Bakalářská práce v elektronické podobě