F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal hemzal@physics.muni.cz MKP — příklad sestavení v ID, Poissonova rovnice Uvažujme Poissonovu rovnici L(u) = Au - 4np = 0. Je-li u = ^2iU[i]N[i], dostáváme pro Galérkinovy podmínky 3- / N\j]Ll^u[i]N[i]ján = ^u[i] J N[j]V(VN[i\)dn-A7T J N\j]páíl = 0, kde ovšem můžeme psát N\j]V(VN[i]) = V(N[j] V N [i]) — VAT [7] VAT [i] a na první člen použít Gaussovu větu. Dostaneme j : ľ«[i] í N[j]VN[i\dW ~y2u[i\ í VA^[j]VA^[ž]d^ = 4tt í N[j]pdn, i JdťL i Jn Jn což je algebraická soustava Ku = b, kde příspěvek objemových členů je Kíj= í VN[j]VN[i\díl bj = An í N[j]pdíl. Omezíme-li se Dirichletovými a homogenními Neumannovými podmínkami, nepřináší okrajový integrál žádný vklad. Přesněji, triviální Neumannova podmínka nevyžaduje žádnou modifikaci rovnic, podmínka Dirichletova svými hodnotami vstoupí do integrálu u okrajových elementů, ale jelikož je hodnota na hranici zadána, všechny takové členy se převádí jako konstanta na pravou stranu (v Dirichletovsky zadaných uzlech rovnice nesestavujeme). F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal hemzal@physics.muni.cz V uvedeném zápisu by sestavení rovnic procházelo ve velké smyčce všechny uzly, v každém cyklu by se dohledávaly všechny elementy, které uzel obsahují. To není praktické, zejména ve vyšších dimenzích výstup z generátoru sítě obsahuje primárně seznam elementů a jeho neustálé prohledávání pro aktuální uzel by bylo extrémně zdlouhavé. Budeme proto matici soustavy a vektor pravé strany sestavovat ve smyčce přes jednotlivé elementy: v ID máme úsečky T(x[i], x[j]) a na nich jednotlivé větve aproximačních funkcí, NT[j]= X~X[%}^ NT[il= X~X[J X[J\~X[l\ X[l\~X[J\ Na každém elementu budou existovat dva typy příspěvků: dva diagonální, typu f 2 1 fxW 1 K33 += / (VA/r[j]) da: = ——-—2 / da: = —--— > 0 pro x[j] > x[i JT {X[JÍ -X[l\y Jx[{] X[j\-X[l\ a mimodiagonální, typu rx[j] i rxiJ} i Kíj += / VNT[j]VNT[i\dx = --—-—2 / da: =----- < 0 pro x[j] > x[i]. {X[j\ -X[l\y JXU] x[j\-x[i\ ' X\l\ Základní kontrolou našeho sestavení je symetrie matice K, v obecnějším případě bývá matice soustavy před zavedením okrajových podmínek hermiteovská. V praxi pořadí uzlů neřešíme a ve jmenovateli použijeme délku elementu It = \x[j] — x[i]\. F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hcmzal hemzal@physics.muni.c z Komponenty vektoru pravé strany b j = 47ľ / N[j]pdil se rovněž rozpadají na příspěvky v rámci jednotlivých elementů. V souladu s myšlenkou MKP můžeme uvažovat p = p[i]N[i], potom 6j=4tt^^[z] í N\j]N[i\d x dá na elementu T(I, J) příspěvky b[J] += 4tt 6 3 lT b[I] += 47T p[J] p\£ 6 3 Vypišme si všechny získané koeficienty ve zjednodušením případě ekvidistantního kroku s: 1 2 1 Kj~lj = ~- K33 = ~ Kj+lj = " c t 8 Pb'-l]+4p[y]+pb" + l]' 6 S e z pravé strany rozeznáváme v sestavené matici diskretizovaný Laplasián, na pravé straně zbývá váhovaný průměr hodnot v okolí aktuálního uzlu, přičemž většina váhy spočívá na hodnotě v aktuálním uzlu. Důvod této změny (ve FD by vystupoval pouze aktuální uzel) tkví v Galerkinově optimalizaci - rovnici se snažíme plnit nejen v uzlových bodech, ale i v oblasti mezi nimi. 4tt F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal hemzal@physics.muni.cz Celé sestavení FEM v našem případě vypadá takto: vytvoření prázdné matice K o velikosti počtu vnitřních uzlů vytvoření prázdného vektoru b o délce počtu vnitřních uzlů pro každý element T(7, J) : It = \x[I] — x[J] | pro vnitřní uzly I, J elementu T: la) Kn += l/lT 2a) 6/ += 47T 6^3 lb) Kjj += l/lT 2b) 6j+=4tt 6^3 lc) Ku += -1/h ld) Kn += -1//T