Chalmers


Datorstudioövning 1 - 2: Assemblera hus

På dagens datorstudioövning är uppgiften att med din egen assemblerare från DTA beräkna matriser och vektorer för "husproblemet", samt jämföra med resultaten du beräknat för hand enligt betinget. Detta är en god kontroll, dels av att du förstår tekniken med assemblering, dels av att dina rutiner fungerar som de skall innan vi går vidare.

Vi skall inte rita och triangulera "huset" i pdetool, utan istället för hand skapa matriserna p (point matrix), t (triangle matrix), och e (edge matrix). De två första matriserna har du nog en viss känsla för hur de är konstruerade, medan den tredje kanske inte är lika bekant. Så att du får göra detta för hand är både avsett som en repetition av datastrukturerna i p och t, samt att lära dig lite mer om e.

Vi börjar med p: denna matris har storleken 2 x nnodes (där nnodes betecknar antalet noder i trianguleringen). Kolonn j innehåller x- och y-koordinaterna (i rad 1 respektive rad 2) för nod nummer j i trianguleringen. Enligt figuren ovan har vi sju stycken noder, och om vi (som i figuren) väljer att lägga origo i husets nedre vänstra hörn kan vi skapa p enligt:


  >> p = zeros(2, 7);    % initiera p
  >> p(:, 1) = [0; 0]    % nod 1: x = 0; y = 0
  >> p(:, 2) = [1; 0]    % nod 2: x = 1; y = 0

Mata nu själv in data för övriga fem noder i p.

Vi fortsätter med t: denna matris har storleken 4 x ntri (där ntri betecknar antalet trianglar i trianguleringen). Raderna 1-3 i kolonn j innehåller nodnumren för de tre hörnen i triangel nummer j, genomlöpta i moturs följd. (Däremot spelar det inte någon roll vilket av hörnen som kommer "först", d.v.s., i rad 1.) Rad 4 anger i vilket delområde triangeln ligger. I detta exempel tänker vi oss att vi endast har ett delområde, hela Omega, vilket har delområdesnummer 1. Notera att området utanför Omega per definition antas ha delområdesnummer 0. (Kommentar: Som exempel på ett fall med flera delområden kan vi tänka oss att "övervåningen", d.v.s. trianglarna 4 och 5, utgör ett delområde (nr. 1, säg), och "nedervåningen", d.v.s. trianglarna 1, 2, 3, och 6, utgör ett annat delområde (nr. 2). Vi skulle då på ett enkelt sätt kunna definiera olika uttryck för t.ex. värmekällan f beroende på om en given triangel hör till delområde 1 eller 2.) Enligt figuren ovan har vi sex stycken trianglar, och vi kan skapa t enligt:


  >> t = zeros(4, 6);        % initiera t
  >> t(:, 1) = [1; 6; 7; 1]  % triangel 1; hörnnoder: 1, 6, 7;
                             %             delområdesnummer: 1
  >> t(:, 2) = [6; 2; 7; 1]  % triangel 2; hörnnoder: 2, 6, 7;
                             %             delområdesnummer: 1

Mata nu själv in data för övriga fyra trianglar i t.

Vi skapar slutligen e: denna matris har storleken 7 x nedges. Här betecknar nedges antalet triangelsidor som ligger på ett randsegment. Randsegment är kurvor som skiljer två delområden åt. I vårt fall har vi endast ett (inre) delområde, hela Omega, och (unionen av) randsegmenten sammanfaller därför med randen till Omega (som skiljer Omega, delområde 1, från området utanför Omega, "delområde" 0). (Kommentar: Om vi istället, enligt kommentaren ovan, tänker oss att "övervåningen" och "nedervåningen" utgör två olika delområden skulle vi få ett extra randsegment: nämligen linjen som sammanbinder noderna 3 och 5. I detta fall skulle e-matrisen alltså ha två extra kolonner, motsvarande de två triangelsidor som ligger på detta, inre, randsegment. Detta randsegment är ett exempel på en så kallad inre rand.) Enligt figuren tänker vi oss att vi har fem stycken randsegment: randsegment nummer 1 är "golv"; randsegment nummer 2 och 5 är "väggar"; randsegment nummer 3 och 4 utgör "taket". Eftersom det ligger två stycken triangelsidor på randsegment nummer 1, och endast en triangelsida på vart och ett av de övriga fyra randsegmenten, blir nedges = 2 + 1 + 1 + 1 + 1 = 6, varför e blir en 7 x 6-matris. Notera också (se figur) att vi har infört en orientering av varje randsegment. Den är vald så att området, Omega, hela tiden ligger till vänster om dig då du "promenerar" längs randen i orienteringens riktning.

Vi inför nu (godtyckligt) följande numrering av triangelsidorna som ligger på ett randsegment:


  1. sidan som sammanbinder nod 1 och nod 6
  2.          "                 6    "    2
  3.          "                 2    "    3
  4.          "                 3    "    4
  5.          "                 4    "    5
  6.          "                 5    "    1

För var och en av dessa sex sidor skall motsvarande kolonn i e innehålla följande information:

  rad 1: nodnummer för sidans "start-nod"
  rad 2: nodnummer för sidans "slut-nod"

    (Vilken som är "start-nod" och vilken som är "slut-nod"
     hänger förstås samman med orienteringen hos det rand-
     segment varpå sidan ligger.)

  rad 3: parametervärde i sidans "start-nod"
  rad 4: parametervärde i sidans "slut-nod"

    (Dessa använder vi oss inte av men om du är nyfiken
     så fungerar det så här: varje randsegment antas vara
     parametriserat med en parameter som är proportionell
     mot båglängden på ett sådant sätt att den "första" noden
     på randsegmentet har parametervärdet 0 och den "sista" har
     parametervärdet 1. När du rör dig längs randsegmentet i
     orienteringens riktning växer alltså parametervärdet från
     0 till 1: t.ex. är det 0.5 när du gått halvvägs!)

  rad 5: numret på det randsegment varpå sidan ligger

  rad 6: numret på delområdet till vänster om sidan
  rad 7: numret på delområdet till höger om sidan

     (När du rör dig i orienteringens riktning.)

Vi kan därmed skapa e enligt:

  >> e = zeros(7, 6);                   % initiera e
  >> e(:, 1) = [1; 6; 0; 0.5; 1; 1; 0]  % triangelsida 1
  >> e(:, 2) = [6; 2; 0.5; 1; 1; 1; 0]  % triangelsida 2

Mata nu själv in data för övriga fyra triangelsidor i e. Notera att delområdet till vänster alltid har nummer 1 medan "delområdet" till höger, utanför Omega, har nummer 0. Rad 3-4 kan du i princip sätta till vad som helst: de används som sagt ej i våra program!

Nu har du skapat matriserna p, t, och e, vilka beskriver trianguleringen av Omega, och vilka behövs som indata till din assemblerare. Men innan du anropar denna måste du även specificera data till problemet i de funktionsfiler som assembleraren i sin tur anropar. I husproblemet så är diffusionskoefficienten och källtermen konstanta (=1) i hela Omega, men tänk på att vi har olika randdata på "golv" (randsegment 1) och "väggar och tak" (randsegment 2-5). Var mycket noggrann när du gör detta, och glöm aldrig bort att kontrollera samtliga "data-filer" innan du startar en ny beräkning.

Anropa nu din(a) assemblerare (obs: inte din "solver"; vi skall inte lösa differentialekvationen, endast assemblera matriser och vektorer), och beräkna 7x7-matriserna M (massmatris), A (styvhetsmatris) och K ("randbidragsmatris"), samt 7x1-vektorerna F (lastvektor) och G ("randbidragsvektor"). Jämför med resultaten från din handräkning!

A, F och G bör stämma överens exakt.

Om du i din kod för att beräkna dubbelintegralerna i massmatriselementen använt kvadratur baserad på integrandens värde i triangelsidornas mittpunkter bör också M bli exakt, men om du använt nodkvadratur får du istället den lumpade massmatrisen; prova i detta fall att för varje rad i din handräknade massmatris addera samtliga element i denna rad, och kontrollera därefter att summan stämmer överens med diagonalelementet i motsvarande rad i den lumpade massmatrisen. Den exakta och den lumpade massmatrisen ser ut som följer:


           /                                     \
           |                                     |
           | 1/24    0    0    0 1/96 1/96 1/48  |
           |      1/24 1/96    0    0 1/96 1/48  |
           |           1/24 1/96    0    0 1/48  |
       M = |                1/24 1/96    0 1/48  |
           |    symm.            1/24    0 1/48  |
           |                          1/24 1/48  |
           |                                1/8  |
           |                                     |
           \                                     /


           /                                     \
           |                                     |
           | 1/12    0    0    0    0    0    0  |
           |      1/12    0    0    0    0    0  |
           |           1/12    0    0    0    0  |
  M_lump = |                1/12    0    0    0  |
           |    symm.            1/12    0    0  |
           |                          1/12    0  |
           |                                1/4  |
           |                                     |
           \                                     /

För att få ut K exakt från din assemblerare krävs att du använder antingen Simpson's formel eller Gauss' 2-punktsformel vid beräkning av kurvintegralerna över randen till Omega. Om du inte modifierat denna del av programmet i den kod du fått utgå ifrån, är det troligt att trapetsregeln används och i dessa fall blir alltså K ej exakt. Titta efter i din assemblerare vilken kvadraturregel som är "i bruk", och jämför sedan med relevant variant nedan:

     /                                                                   \
     |                                                                   |
     | 1/6   0                   0             0                1/12 0 0 |
     |     1/6                1/12             0                   0 0 0 |
     |         (1/6)+1/(3*sqrt(2)) 1/(6*sqrt(2))                   0 0 0 |
 K = |                             2/(3*sqrt(2))       1/(6*sqrt(2)) 0 0 |
     |    symm.                                  (1/6)+1/(3*sqrt(2)) 0 0 |
     |                                                               0 0 |
     |                                                                 0 |
     |                                                                   |
     \                                                                   /


          /                                                               \
          |                                                               |
          | 1/4   0                   0         0                   0 0 0 |
          |     1/4                   0         0                   0 0 0 |
          |         (1/4)+1/(2*sqrt(2))         0                   0 0 0 |
 K_trap = |                             1/sqrt(2)                   0 0 0 |
          |    symm.                              (1/4)+1/(2*sqrt(2)) 0 0 |
          |                                                           0 0 |
          |                                                             0 |
          |                                                               |
          \                                                               /


         /                                                                   \
         |                                                                   |
         | 1/8   0                   0             0                 1/8 0 0 |
         |     1/8                 1/8             0                   0 0 0 |
         |         (1/8)+1/(4*sqrt(2)) 1/(4*sqrt(2))                   0 0 0 |
 K_mid = |                             1/(2*sqrt(2))       1/(4*sqrt(2)) 0 0 |
         |    symm.                                  (1/8)+1/(4*sqrt(2)) 0 0 |
         |                                                               0 0 |
         |                                                                 0 |
         |                                                                   |
         \                                                                   /

Extra tester du kan göra:

1. Litet hus. Prova att förminska huset med en faktor två (i längdskala) genom kommandot >> p = p/2! Räkna ut nya matriser och vektorer för det lilla huset (i Matlab). Jämför nu den gamla och den nya massmatrisen. Hur skiljer de sig åt? Förklara! Gör samma sak med styvhetsmatrisen. Kommentar?

2. Radsummor i M. Beräkna kolonnvektorn med radsummor av elementen i massmatrisen, d.v.s., kolonnvektorn där elementet i rad i är summan av samtliga element i motsvarande rad i M. (Jfr. lumpad massmatris.) Tips: Använd kommandot >> radsummor = sum(M, 2). Jämför med lastvektorn F. Kan du förklara varför dessa är lika i detta fall (med f=1)? Ge en geometrisk tolkning av elementen i radsummor som volymer. Ledtråd: Summan av alla tältfunktioner är den funktion som är identiskt lika med 1.

3. Arean av Omega. Räkna nu ut summan av samtliga element i kolonnvektorn från uppgift 2. Tips: Använd kommandot >> sum(radsummor, 1). Jämför resultatet med arean av Omega. Förklara!

4. Radsummor i A. Beräkna kolonnvektorn med radsummor av elementen i styvhetsmatrisen, analogt med uppgift 2. Vad blir resultatet? Varför?

Last modified: Tue Sep 10 10:53:51 MET DST 2002