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