function [S, M, R, v, r] = PoissonAssembler2D(p,e,t,h);

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

S = zeros(n,n);     % Initiera styvhetsmatris.
M = zeros(n,n);     % Initiera massmatris.
R = zeros(n,n);     % Initiera randmatris.
v = zeros(n,1);     % Initiera lastvektor.
r = zeros(n,1);     % Initiera randvektor.

for el = 1:ntri

  nodes = t(1:3,el);
  coords = p(:,nodes);
  
  Me = ElementMassmatrix(h);
  Se = ElementStiffnessmatrix(h);
  ve = ElementLoadvector(coords,h);
 
  M(nodes,nodes) = M(nodes,nodes) + Me;
  S(nodes,nodes) = S(nodes,nodes) + Se;
  v(nodes) = v(nodes) + ve; 

end

%
% Bidrag från randen, OBS : ÄNDRA INTE
%

for bel = 1:size(e,2)

    nodes = e(1:2,bel);
    coords = p(:,nodes);
    g_N = 0.0; 
    g_D = 0.0;
    gamma = 1e5; 
    sidelength = norm(coords(:,1)-coords(:,2));
    phi = [.5 .5];
    R(nodes,nodes) = R(nodes,nodes) + gamma*phi'*phi*sidelength; 
    r(nodes) = r(nodes) + (gamma*g_D - g_N)*phi'*sidelength;

end

% subroutines -----------------------------------------------------------------

function Me = ElementMassmatrix(h)
%
% Komplettera med de rätta värdena på elementmassmatrisen Me.
%
Me = [ 0.0 , 0.0 , 0.0;
       0.0 , 0.0 , 0.0;
       0.0 , 0.0 , 0.0 ];
%
function Se = ElementStiffnessmatrix(h)
%
% Komplettera med de rätta värdena på elementstyvhetsmatrisen Se.
%
Se = [ 0.0 , 0.0 , 0.0;
       0.0 , 0.0 , 0.0;
       0.0 , 0.0 , 0.0 ];
%
function ve = ElementLoadvector(coords,h)
%
% Använd kvadratur för att beräkna elementlastvektorn ve.
%
triangelarea = 0.0;
%
x = coords(1,:); 
y = coords(2,:);
ve = [ 0.0 ; 0.0 ; 0.0 ];
%
% Lasten f.
%
function load = f(x,y)
load = 0.0;