K2/Bt2 - TMK: Tillämpad matematik K (TMA 682) - 2003/2004

Datorstudioövning 5: Kvadratur i 2D. Styvhetsmatris. Robinrandvillkor.



Dagens uppgift är att bygga ut programmen PoissonSolver2D.m och PoissonAssembler2D.n till fullständiga FEM-lösare för stationära reaktions-diffusonsproblem med Dirichlet villkor (Robinrandvillkor) i två rumsdimensioner:

-&Delta u+u=f, i D=(0,1)x(0,1),
u=0, på randen av D.

Börja med att ladda ner filerna PoissonSolver2D.m och PoissonAssembler2D.m till ditt hembibliotek. Starta sedan Matlab.

Dagens uppgift är att bygga ut programmet PoissonSolver2D.m till en komplett FEM-lösare för stationära reaktions-diffusionsproblem med homogena Dirichletvillkor i två rumsdimensioner.

Börja med att variationsformulera och diskretisera så att du får fram det relevanta ekvationssystemet på matrisform. Kom ihåg att huvudmomenten för detta är: Multiplicera med en testfunktion och integrera över domänen. Partialintegera och ansätt en lösningsapproximation i form av en summa av tältfunktioner och koefficienter. Välj testfunktioner som tältfunktionerna och skriv slutligen på matrisform.

Huvudmomenten i PoissonSolver.m är: 

  1. Triangulering av domänen [0,1]x[0,1].

  2. Assemblering av matriser (S, M, R) och vektorer (v, r).

  3. Lösning av det linjära ekvationssystemet U = A\b och visualisering.
Dessa steg är implementerade på följande sätt:

1. Först skapas en triangulering i huvudprogrammet PoissonSolver2D.m mha. Matlabs meshgenerator poimesh (p,e,t). Här är p, e och t matriser, som beskriver trianguleringens nodkoordinater, konnektivitet och kanter.

2. Därefter assembleras mass- och styvhetsmatriserna i PoissonAssembler2D.m. Även lastvektorn assembleras här. För att hantera olika typer av randvillkor, sk. Robin randvillkor finns även en randmatris R och en randvektor r. Du behöver dock inte göra något åt dessa eftersom de leder till ett homogent randvillkor av Dirichlet typ i vår kod.

3. Lösning av det linjära ekvationssystemet och plottning av lösningen görs slutligen i huvudprogrammet.

I likhet med förra studioövningen måste du komplettera koden för att få en fungerande FEM-lösare. Du skall skriva in elementmatriserna Me och Se för massmatrisen och styvhetsmatrisen, samt elementlastvektorn ve. Använd resultat och metod från förra veckan. För lastvektorns del skall du använda nodkvadratur för att approximativt beräkna integralen. Se kapitel 76.7 formel (76.14). Börja med en konstant last, tex.
f(x, y) = 1.0, vilket gör att du lätt kan kontrollera om ve blir rätt (det borde bli vektorn ve = h*[h,h,h]/6). Prova sedan med en last som har tydliga karaktäristika i x och y-led, tex. f(x, y) = exp(y) sin 7x. Fundera på hur detta avspeglas i lösningen på skärmen.

Lycka till!

Extra uppgift:

Läs igenom materialet om Robinrandvillkor och utöka ditt program så att det kan hantera olika randvillkor.

Tillbaka till kurshemsidan.


Editor: Fredrik Bengzon
Last modified: 2003-09-26