function [A,M] = matrix_assembly(Elemente,Knoten)
% Diese Funktion berechnet für die DGL -d*nabla(u) + cu = f im Inneren, 
% u = 0 auf dem Rand die Steifigkeitsmatrix A und der Massematrix M
%
% Input:  Elemente aus txt file
%         Knoten aus txt file
%
% Output: A Matrix (Steifigkeitsmatrix)
%         M Matrix (Massematrix)


[n1,~] = size(Elemente);
[n2,~] = size(Knoten);
A = sparse(n2,n2);            % Steifigkeitsmatrix
M = sparse(n2,n2);            % Massematrix

% Berechnet in Aufgabe 23 Blatt 9 für Referenzdreieck
M0 = [1/12,1/24,1/24;1/24,1/12,1/24;1/24,1/24,1/12];
Axx = 0.5*[1,-1,0;-1,1,0;0,0,0];
Axy = 0.5*[1,0,-1;-1,0,1;0,0,0];
Ayx = 0.5*[1,-1,0;0,0,0;-1,1,0];
Ayy = 0.5*[1,0,-1;0,0,0;-1,0,1];

for j=1:n1
    
    % Transformation vom Einheitsdreieck auf das j-te Dreieck der Triangulierung
    L = [Knoten(Elemente(j,2),1) - Knoten(Elemente(j,1),1),Knoten(Elemente(j,3),1) - Knoten(Elemente(j,1),1);
        Knoten(Elemente(j,2),2) - Knoten(Elemente(j,1),2),Knoten(Elemente(j,3),2) - Knoten(Elemente(j,1),2)];
    C = inv(L);
    Determinante = abs(det(L));   % Funktionaldeterminante
    % Ende Transformation
    
    % Berechnung der lokalen Matrizen, siehe Aufgabe 26
    Mloc = Determinante*M0;
    Aloc = Determinante*((C(1,1)^2+C(1,2)^2)*Axx + (C(2,1)^2+C(2,2)^2)*Ayy + (C(1,1)*C(2,1)+C(1,2)*C(2,2))*(Axy+Ayx));
    
    % Beschreibung der globalen Matrizen A und M
    A(Elemente(j,1:3),Elemente(j,1:3)) = A(Elemente(j,1:3),Elemente(j,1:3)) + Aloc;
    M(Elemente(j,1:3),Elemente(j,1:3)) = M(Elemente(j,1:3),Elemente(j,1:3)) + Mloc;
    
end
end