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.
|