Kf1 - Differentialekvationer och tekniska beräkningar del A, TMA 225 - 2003

Datorstudioövning 7: Kvadratur (2D). Massmatris.

Nedan finns koden till en funktion som med nodkvadratur beräknar en approximation till massmatrisen. Inargument är de två matriserna p och t som beskriver din triangulering. Se ``Problems Week 5'' om du inte minns hur dessa ser ut. Skapa själv en funktion som heter Min2DMassmatrisAssemblerare.m där du skriver in koden. Tänk noga igenom varje steg så att du förstår hur koden fungerar.


  function M = Min2DMassmatrisAssemblerare(p, t);

  nnodes = size(p,2);  % Antal noder. (= antal kolonner i p)
  ntri = size(t,2);    % Antal trianglar. (= antal kolonner i t)

  M = zeros(nnodes, nnodes);  % Initiera massmatris.

  for n = 1:ntri

    % Ta reda på hörnens nodnummer:  
  
    i = t(1,n);  % Nodnummer för hörn nummer 1 i triangel nummer n.
    j = t(2,n);  % Nodnummer för hörn nummer 2 i triangel nummer n.
    k = t(3,n);  % Nodnummer för hörn nummer 3 i triangel nummer n.
                 %
                 % 1 -> 2 -> 3 moturs riktning.
  
    % Ta reda på hörnens koordinater:                      
                       
    xnodi = p(1,i);  % x-koordinat för nod nummer i.
    xnodj = p(1,j);  % x-koordinat för nod nummer j.
    xnodk = p(1,k);  % x-koordinat för nod nummer k.
  
    ynodi = p(2,i);  % y-koordinat för nod nummer i.
    ynodj = p(2,j);  % y-koordinat för nod nummer j.
    ynodk = p(2,k);  % y-koordinat för nod nummer k.
                       
    % Räkna ut triangelns area:
    % Jämför Problem 4 i ``Problems Week 6''.
  
    triangelarea = ((xnodj - xnodi)*(ynodk - ynodi) - ...
                    (xnodk - xnodi)*(ynodj - ynodi))/2;

    % ... gör att man kan fortsätta skriva på nästa rad.

    % Räkna ut värdena på phi_i, phi_j, och phi_k, i triangelns hörn:
  
    phi_inodi = 1;       % Värdet i nod i av phi_i.
    phi_inodj = 0;       % Värdet i nod j av phi_i.
    phi_inodk = 0;       % Värdet i nod k av phi_i.
  
    phi_jnodi = 0;       % Värdet i nod i av phi_j.
    phi_jnodj = 1;       % Värdet i nod j av phi_j.
    phi_jnodk = 0;       % Värdet i nod k av phi_j.
  
    phi_knodi = 0;       % Värdet i nod i av phi_k.
    phi_knodj = 0;       % Värdet i nod j av phi_k.
    phi_knodk = 1;       % Värdet i nod k av phi_k.

    % Bidragen till massmatrisen från triangel nummer n
    % beräknas med hörnkvadratur:

    M(i,i) = M(i,i) + ...
        (phi_inodi*phi_inodi + ...
         phi_inodj*phi_inodj + ...
         phi_inodk*phi_inodk)*triangelarea/3;
   
    M(i,j) = M(i,j) + ...
        (phi_jnodi*phi_inodi + ...
         phi_jnodj*phi_inodj + ...
         phi_jnodk*phi_inodk)*triangelarea/3;
   
    M(i,k) = M(i,k) + ...
        (phi_knodi*phi_inodi + ...
         phi_knodj*phi_inodj + ...
         phi_knodk*phi_inodk)*triangelarea/3;
   
    M(j,i) = M(j,i) + ...
        (phi_inodi*phi_jnodi + ...
         phi_inodj*phi_jnodj + ...
         phi_inodk*phi_jnodk)*triangelarea/3;
   
    M(j,j) = M(j,j) + ...
        (phi_jnodi*phi_jnodi + ...
         phi_jnodj*phi_jnodj + ...
         phi_jnodk*phi_jnodk)*triangelarea/3;
   
    M(j,k) = M(j,k) + ...
        (phi_knodi*phi_jnodi + ...
         phi_knodj*phi_jnodj + ...
         phi_knodk*phi_jnodk)*triangelarea/3;
   
    M(k,i) = M(k,i) + ...
        (phi_inodi*phi_knodi + ...
         phi_inodj*phi_knodj + ...
         phi_inodk*phi_knodk)*triangelarea/3;
   
    M(k,j) = M(k,j) + ...
        (phi_jnodi*phi_knodi + ...
         phi_jnodj*phi_knodj + ...
         phi_jnodk*phi_knodk)*triangelarea/3;
   
    M(k,k) = M(k,k) + ...
        (phi_knodi*phi_knodi + ...
         phi_knodj*phi_knodj + ...
         phi_knodk*phi_knodk)*triangelarea/3;

  end

Skapa nu matriserna p och t som beskriver trianguleringen i Problem 5 i ``Problems Week 6'':

  >> p(:,1) = [0; 0];
  >> ...
  >> p(:,5) = ...
  >> t(:,1) = [1; 2; 5];
  >> ...
  >> t(:,3) = ...

och beräkna (den approximativa) massmatrisen:

  >> M = Min2DMassmatrisAssemblerare(p, t)

Man kan visa (se Problem 5 (c*) i ``Problems Week 6'') att den approximation som fås med nodkvadratur inte är något annat än den ``lumpade'' massmatrisen! Jämför därför med ditt resultat från Problem 5 (b) i ``Problems Week 6''. Uppgiften är nu att byta kvadraturformel i Min2DMassmatrisAssemblerare från nodkvadratur till kvadratur baserad på integrandens värde i mittpunkterna på triangelsidorna. Denna kvadraturformel är exakt för integrander som är polynom av grad 2 och ger alltså exakt resultat vid beräkning av massmatrisen! Prova nu därför att beräkna:

  >> M = MinNya2DMassmatrisAssemblerare(p, t)

och jämför med resultatet från Problem 5 (a) i ``Problems Week 6''.

Lycka till!

Tillbaka till kurshemsidan.


Editor: Georgios Foufas
Last modified: 2003-04-24