function [A, Mc, R, b, rv] = MyFirst2DPoissonAssembler(p, e, t); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% % %% Syfte: Assemblerar matriser och vektorer för att, med cG(1)-metoden, lösa % %% differentialekvationen % %% % %% -div(a(x,y)*grad(u(x,y))) + c(x,y)*u(x,y) = f(x,y), i Omega, % %% % %% med Robinrandvillkor % %% % %% -a*Du/Dn = gamma*(u - g_D) + g_N, på randen, Gamma, till Omega. % %% % %% % %% Indata: % %% % %% p: "Point matrix" (Matris med dimension 2 x nnodes) % %% e: "Edge matrix" (Matris med dimension 7 x nedges) % %% t: "Triangle matrix" (Matris med dimension 4 x ntri) % %% % %% % %% Utdata: % %% A: styvhetsmatris (Matris med dimension nnodes x nnodes) % %% Mc: massmatris (Matris med dimension nnodes x nnodes) % %% R: randmatris (Matris med dimension nnodes x nnodes) % %% % %% b: lastvektor (Vektor med dimension nnodes x 1) % %% rv: randvektor (Vektor med dimension nnodes x 1) % %% % %% % %% % %% Anropade funktioner: a_2D.m, c_2D.m, f_2D.m, gamma_2D.m, g_D_2D.m, och g_N_2D.m % %% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Här börjar MyFirst2DPoissonAssembler. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nnodes = size(p,2); % Antal noder. (= antal kolonner i p) nedges = size(e,2); % Antal triangelsidor som ligger på ett randsegment. % (= antal kolonner i e) ntri = size(t,2); % Antal trianglar. (= antal kolonner i t) A = zeros(nnodes, nnodes); % Initiera styvhetsmatris. Mc = zeros(nnodes, nnodes); % Initiera massmatris. R = zeros(nnodes, nnodes); % Initiera randmatris. b = zeros(nnodes, 1); % Initiera lastvektor. rv = zeros(nnodes, 1); % Initiera randvektor. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Vi adderar först bidrag fr.o.m. triangel 1 t.o.m. triangel ntri. % % % % Triangel nummer n, vars hörn har nodnummer: % % % % i = t(1,n) % % j = t(2,n) % % k = t(3,n) % % % % bidrar till följande element i 'A' och 'Mc': % % % % A(i,i) A(i,j) A(i,k) Mc(i,i) Mc(i,j) Mc(i,k) % % A(j,i) A(j,j) A(j,k) Mc(j,i) Mc(j,j) Mc(j,k) % % A(k,i) A(k,j) A(k,k) Mc(k,i) Mc(k,j) Mc(k,k) % % % % och till följande element i 'b': % % % % b(i) % % b(j) % % b(k) % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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. % Beräkna koordinaterna för triangelsidornas mittpunkter: xmittij = (xnodi + xnodj)/2; % Mittpunkt på sida ij xmittik = (xnodi + xnodk)/2; % Mittpunkt på sida ik xmittjk = (xnodj + xnodk)/2; % Mittpunkt på sida jk ymittij = (ynodi + ynodj)/2; ymittik = (ynodi + ynodk)/2; ymittjk = (ynodj + ynodk)/2; % Räkna ut triangelns area: triangelarea = ((xnodj - xnodi)*(ynodk - ynodi) - ... (xnodk - xnodi)*(ynodj - ynodi))/2; % ... gör att man kan fortsätta skriva på nästa rad. % I vilket delområde ligger triangel nummer n? delomradesnummer = t(4,n); % Rad 4 i t innehåller numret på delområdet % i vilket triangel nummer n ligger. % Räkna ut diffusionskoefficienten och källtermen i triangelns hörn, % samt hastighetskonstanten i triangelsidornas mittpunkter: % (Vi skickar med numret på det delområde i vilket triangeln ligger % som en parameter, för att på ett enkelt sätt kunna specificera olika % uttryck för a, c, och f i skilda delområden.) anodi = a_2D(xnodi,ynodi,delomradesnummer); % Värdet i nod i av a. anodj = a_2D(xnodj,ynodj,delomradesnummer); % Värdet i nod j av a. anodk = a_2D(xnodk,ynodk,delomradesnummer); % Värdet i nod k av a. fnodi = f_2D(xnodi,ynodi,delomradesnummer); % Värdet i nod i av f. fnodj = f_2D(xnodj,ynodj,delomradesnummer); % Värdet i nod j av f. fnodk = f_2D(xnodk,ynodk,delomradesnummer); % Värdet i nod k av f. cmittij = c_2D(xmittij,ymittij,delomradesnummer); % Värdet av c mitt på sida ij. cmittik = c_2D(xmittik,ymittik,delomradesnummer); % Värdet av c mitt på sida ik. cmittjk = c_2D(xmittjk,ymittjk,delomradesnummer); % Värdet av c mitt på sida jk. % 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. % Beräkna värdena på phi_i, phi_j, och phi_k, mitt på triangelsidorna: phi_imittij = 1/2; % Värdet mitt på sida ij av phi_i. phi_imittik = 1/2; % Värdet mitt på sida ik av phi_i. phi_imittjk = 0; % Värdet mitt på sida jk av phi_i. phi_jmittij = 1/2; % Värdet mitt på sida ij av phi_j. phi_jmittik = 0; % Värdet mitt på sida ik av phi_j. phi_jmittjk = 1/2; % Värdet mitt på sida jk av phi_j. phi_kmittij = 0; % Värdet mitt på sida ij av phi_k. phi_kmittik = 1/2; % Värdet mitt på sida ik av phi_k. phi_kmittjk = 1/2; % Värdet mitt på sida jk av phi_k. % Räkna ut (de konstanta) gradienterna av phi_i, phi_j, och phi_k: Dphi_iDx = (ynodj - ynodk)/(2*triangelarea); Dphi_iDy = (xnodk - xnodj)/(2*triangelarea); gradphi_i = [Dphi_iDx; Dphi_iDy]; Dphi_jDx = (ynodk - ynodi)/(2*triangelarea); Dphi_jDy = (xnodi - xnodk)/(2*triangelarea); gradphi_j = [Dphi_jDx; Dphi_jDy]; Dphi_kDx = (ynodi - ynodj)/(2*triangelarea); Dphi_kDy = (xnodj - xnodi)/(2*triangelarea); gradphi_k = [Dphi_kDx; Dphi_kDy]; % Hörnkvadratur: % Notera: dot(v,w) beräknar skalärprodukten mellan vektorerna v och w. A(i,i) = A(i,i) + ... % Integrand: a*dot(gradphi_i,gradphi_i) (anodi*dot(gradphi_i,gradphi_i) + ... anodj*dot(gradphi_i,gradphi_i) + ... anodk*dot(gradphi_i,gradphi_i))*triangelarea/3; A(i,j) = A(i,j) + ... % Integrand: a*dot(gradphi_j,gradphi_i) (anodi*dot(gradphi_j,gradphi_i) + ... anodj*dot(gradphi_j,gradphi_i) + ... anodk*dot(gradphi_j,gradphi_i))*triangelarea/3; A(i,k) = A(i,k) + ... (anodi*dot(gradphi_k,gradphi_i) + ... anodj*dot(gradphi_k,gradphi_i) + ... anodk*dot(gradphi_k,gradphi_i))*triangelarea/3; A(j,i) = A(j,i) + ... (anodi*dot(gradphi_i,gradphi_j) + ... anodj*dot(gradphi_i,gradphi_j) + ... anodk*dot(gradphi_i,gradphi_j))*triangelarea/3; A(j,j) = A(j,j) + ... (anodi*dot(gradphi_j,gradphi_j) + ... anodj*dot(gradphi_j,gradphi_j) + ... anodk*dot(gradphi_j,gradphi_j))*triangelarea/3; A(j,k) = A(j,k) + ... (anodi*dot(gradphi_k,gradphi_j) + ... anodj*dot(gradphi_k,gradphi_j) + ... anodk*dot(gradphi_k,gradphi_j))*triangelarea/3; A(k,i) = A(k,i) + ... (anodi*dot(gradphi_i,gradphi_k) + ... anodj*dot(gradphi_i,gradphi_k) + ... anodk*dot(gradphi_i,gradphi_k))*triangelarea/3; A(k,j) = A(k,j) + ... (anodi*dot(gradphi_j,gradphi_k) + ... anodj*dot(gradphi_j,gradphi_k) + ... anodk*dot(gradphi_j,gradphi_k))*triangelarea/3; A(k,k) = A(k,k) + ... (anodi*dot(gradphi_k,gradphi_k) + ... anodj*dot(gradphi_k,gradphi_k) + ... anodk*dot(gradphi_k,gradphi_k))*triangelarea/3; b(i) = b(i) + ... % Integrand: f*phi_i (fnodi*phi_inodi + fnodj*phi_inodj + fnodk*phi_inodk)*triangelarea/3; b(j) = b(j) + ... (fnodi*phi_jnodi + fnodj*phi_jnodj + fnodk*phi_jnodk)*triangelarea/3; b(k) = b(k) + ... (fnodi*phi_knodi + fnodj*phi_knodj + fnodk*phi_knodk)*triangelarea/3; % Mittpunktskvadratur: Mc(i,i) = Mc(i,i) + ... % Integrand: c*phi_i*phi_i (cmittij*phi_imittij*phi_imittij + ... cmittik*phi_imittik*phi_imittik + ... cmittjk*phi_imittjk*phi_imittjk)*triangelarea/3; Mc(i,j) = Mc(i,j) + ... % Integrand: c*phi_j*phi_i (cmittij*phi_jmittij*phi_imittij + ... cmittik*phi_jmittik*phi_imittik + ... cmittjk*phi_jmittjk*phi_imittjk)*triangelarea/3; Mc(i,k) = Mc(i,k) + ... % Integrand: c*phi_k*phi_i (cmittij*phi_kmittij*phi_imittij + ... cmittik*phi_kmittik*phi_imittik + ... cmittjk*phi_kmittjk*phi_imittjk)*triangelarea/3; Mc(j,i) = Mc(j,i) + ... % Integrand: c*phi_i*phi_j (cmittij*phi_imittij*phi_jmittij + ... cmittik*phi_imittik*phi_jmittik + ... cmittjk*phi_imittjk*phi_jmittjk)*triangelarea/3; Mc(j,j) = Mc(j,j) + ... % Integrand: c*phi_j*phi_j (cmittij*phi_jmittij*phi_jmittij + ... cmittik*phi_jmittik*phi_jmittik + ... cmittjk*phi_jmittjk*phi_jmittjk)*triangelarea/3; Mc(j,k) = Mc(j,k) + ... % Integrand: c*phi_k*phi_j (cmittij*phi_kmittij*phi_jmittij + ... cmittik*phi_kmittik*phi_jmittik + ... cmittjk*phi_kmittjk*phi_jmittjk)*triangelarea/3; Mc(k,i) = Mc(k,i) + ... % Integrand: c*phi_i*phi_k (cmittij*phi_imittij*phi_kmittij + ... cmittik*phi_imittik*phi_kmittik + ... cmittjk*phi_imittjk*phi_kmittjk)*triangelarea/3; Mc(k,j) = Mc(k,j) + ... % Integrand: c*phi_j*phi_k (cmittij*phi_jmittij*phi_kmittij + ... cmittik*phi_jmittik*phi_kmittik + ... cmittjk*phi_jmittjk*phi_kmittjk)*triangelarea/3; Mc(k,k) = Mc(k,k) + ... % Integrand: c*phi_k*phi_k (cmittij*phi_kmittij*phi_kmittij + ... cmittik*phi_kmittik*phi_kmittik + ... cmittjk*phi_kmittjk*phi_kmittjk)*triangelarea/3; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Vi adderar nu randbidrag från de triangelsidor som ligger på ett % % (yttre) randsegment. % % % % Triangelsida nummer l, vars "startpunkt" har nodnummer % % % % i = e(1,l) % % % % och vars "slutpunkt" har nodnummer % % % % j = e(2,l) % % % % bidrar (om e(6,l) == 0 eller e(7,l) == 0, d.v.s., om triangelsidan % % ligger på ett yttre randsegment) till följande element % % i 'R': % % % % R(i,i) R(i,j) % % R(j,i) R(j,j) % % % % och till följande element i 'rv': % % % % rv(i) % % rv(j) % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for l = 1:nedges vansterdelomrade = e(6,l); % Numret på delområdet som ligger till % vänster om triangelsida l, när du rör % dig i riktning från nod e(1,l) till % nod e(2,l). hogerdelomrade = e(7,l); if ((vansterdelomrade == 0) | (hogerdelomrade == 0)) % Yttre randsegment! % (Området utanför Omega har delområdesnummer 0.) % Ta reda på ändpunkternas nodnummer: i = e(1,l); % Nodnummer för triangelsida l:s "startpunkt". j = e(2,l); % Nodnummer för triangelsida l:s "slutpunkt". % Ta reda på ändpunkternas koordinater: xnodi = p(1,i); % x-koordinat för nod nummer i ("startpunkt"). xnodj = p(1,j); % x-koordinat för nod nummer j ("slutpunkt"). ynodi = p(2,i); % y-koordinat för nod nummer i ("startpunkt"). ynodj = p(2,j); % y-koordinat för nod nummer j ("slutpunkt"). % Räkna ut triangelsidans längd: sidelength = sqrt((xnodi - xnodj)^2 + (ynodi - ynodj)^2); % På vilket randsegment ligger triangelsida l? randsegmentnummer = e(5,l); % Räkna ut 'gamma', 'g_D', och 'g_N' i triangelsida l:s ändpunkter: % (Vi skickar med numret på det randsegment varpå triangelsidan ligger % som en parameter, för att kunna specificera olika randvillkor på % skilda randsegment.) gammanodi = gamma_2D(xnodi,ynodi,randsegmentnummer); gammanodj = gamma_2D(xnodj,ynodj,randsegmentnummer); g_Dnodi = g_D_2D(xnodi,ynodi,randsegmentnummer); g_Dnodj = g_D_2D(xnodj,ynodj,randsegmentnummer); g_Nnodi = g_N_2D(xnodi,ynodi,randsegmentnummer); g_Nnodj = g_N_2D(xnodj,ynodj,randsegmentnummer); % Räkna ut värdena på phi_i och phi_j i triangelsidans ändpunkter: phi_inodi = 1; phi_inodj = 0; phi_jnodi = 0; phi_jnodj = 1; % Trapetsregeln: R(i,i) = R(i,i) + (gammanodi*phi_inodi*phi_inodi + ... gammanodj*phi_inodj*phi_inodj)/2*sidelength; R(i,j) = R(i,j) + (gammanodi*phi_jnodi*phi_inodi + ... gammanodj*phi_jnodj*phi_inodj)/2*sidelength; R(j,i) = R(j,i) + (gammanodi*phi_inodi*phi_jnodi + ... gammanodj*phi_inodj*phi_jnodj)/2*sidelength; R(j,j) = R(j,j) + (gammanodi*phi_jnodi*phi_jnodi + ... gammanodj*phi_jnodj*phi_jnodj)/2*sidelength; rv(i) = rv(i) + ((gammanodi*g_Dnodi - g_Nnodi)*phi_inodi + ... (gammanodj*g_Dnodj - g_Nnodj)*phi_inodj)/2*sidelength; rv(j) = rv(j) + ((gammanodi*g_Dnodi - g_Nnodi)*phi_jnodi + ... (gammanodj*g_Dnodj - g_Nnodj)*phi_jnodj)/2*sidelength; end % Slut på kontroll om triangelsidan ligger på ett yttre randsegment. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Slut på MyFirst2DPoissonAssembler. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%