function Blatt6

clear all
close all

% Der Code zur Erstellung der Triangulierung, sowie das finden der Dirichlet-Indizes.
addpath('/home/misc/jahn/Lehre/WS1718/Part_WS_17_18/NumPDGL_WS1617/FEM_Code_NPDGL') % hier wird der Pfad angegeben, in welchem das 'distmesh'-Paket liegt
% Polyeder für das Gebiet
pfix = [[-4,-4];[-4,4];[0,4];[0,0];[4,0];[4,-4];[-4,-4]];
% Triangulierung erstellen
[P,T] = distmesh2d(@dpoly,@huniform,0.1,[-5,-5;5,5],pfix,pfix);
% Finde die Dirichlet Knoten
Dirichlet = find(abs(dpoly(P,pfix)) < 1e-3);
% f ist 1 im Teilgebiet D und 0 sonst: verwende f später um zu bestimmen ob ein Dreieck im Teilgebiet ist oder nicht
f = @(p) (p(:,1) >= -3).*(p(:,1) <= -1).*(p(:,2) >= 1).*(p(:,2) <= 3);

% Initialisierung
nt = 1000; % number of timesteps
tau = 1/100; % timestepping
n = length(P); % Anzahl Diskretisierungspunkte
B = spalloc(n,n,5*n); % Matrix B als sparse Matrix
f_val = zeros(3*length(T),3); % Sturktur für das Zwischenspeichern der Lambda_iLambda_j Integrale für die zeitabhängige rechte Seite f
y_part_one = zeros(n,1); % ersten Teil der rechten Seite initialisieren. Volle rechte Seite ist dann 1. Teil + in jedem Zeitschritt der entsprechende 2. Teil

for i = 1:length(T)
    % Wir bestimmen zuerst die Matrix Ainv, womit dann das LGS zur
    % Bestimmung der stk. linearen Funktionen über dem aktuellen Dreieck
    % für verschiedene rechte Seiten gelöst werden kann (bzw. betrachten
    % wir nur Hutfunktionen, sodass die rechte Seite ein Einheitsvektor ist
    % und die Lösung also die entsprechende Spalte von Ainv ist).
    
    % Determinante von A, vgl. Theorie.
    det = P(T(i,2),1)*P(T(i,3),2) + P(T(i,1),1)*P(T(i,2),2) + P(T(i,1),2)*P(T(i,3),1) ... 
        - P(T(i,1),2)*P(T(i,2),1) - P(T(i,3),1)*P(T(i,2),2) - P(T(i,3),2)*P(T(i,1),1);
    % Damit Ainv angeben
    Ainv = 1/det*[P(T(i,2),1)*P(T(i,3),2) - P(T(i,2),2)*P(T(i,3),1), P(T(i,1),2)*P(T(i,3),1) - P(T(i,1),1)*P(T(i,3),2), P(T(i,1),1)*P(T(i,2),2) - P(T(i,1),2)*P(T(i,2),1);...
                  P(T(i,2),2) - P(T(i,3),2), P(T(i,3),2) - P(T(i,1),2), P(T(i,1),2) - P(T(i,2),2);...
                  P(T(i,3),1) - P(T(i,2),1), P(T(i,1),1) - P(T(i,3),1), P(T(i,2),1) - P(T(i,1),1)];
    % Flächeninahlt des Dreiecks - einerseits relevant für die Integrale
    % über Gradient (ist stkw. konstant) - andererseits entspricht 2*vol 
    % der Determinante welche aus der Koordinatentransformation resultiert.
    vol = Inhalt(P(T(i,1),:),P(T(i,2),:),P(T(i,3),:));
              
    % Damit elementweises aufstellen der Matrix B - p1,p2,p3 sind die 3
    % Hutfunktionen und ihre Koeffizienten stehen in der 1.,2.,3. Spalte von Ainv
    % Hier nur der 1. schon bekannte Teil der BLF!
    B(T(i,1),T(i,1)) = B(T(i,1),T(i,1)) + tau*vol*(Ainv(2,1)^2 + Ainv(3,1)^2); % p1 mit p1
    B(T(i,1),T(i,2)) = B(T(i,1),T(i,2)) + tau*vol*(Ainv(2,1)*Ainv(2,2) + Ainv(3,1)*Ainv(3,2)); % p1 mit p2
    B(T(i,1),T(i,3)) = B(T(i,1),T(i,3)) + tau*vol*(Ainv(2,1)*Ainv(2,3) + Ainv(3,1)*Ainv(3,3)); % p1 mit p3
    B(T(i,2),T(i,1)) = B(T(i,2),T(i,1)) + tau*vol*(Ainv(2,2)*Ainv(2,1) + Ainv(3,2)*Ainv(3,1)); % p2 mit p1 - selbe wie p1 mit p2
    B(T(i,2),T(i,2)) = B(T(i,2),T(i,2)) + tau*vol*(Ainv(2,2)^2 + Ainv(3,2)^2); % p2 mit p2
    B(T(i,2),T(i,3)) = B(T(i,2),T(i,3)) + tau*vol*(Ainv(2,2)*Ainv(2,3) + Ainv(3,2)*Ainv(3,3)); % p2 mit p3
    B(T(i,3),T(i,1)) = B(T(i,3),T(i,1)) + tau*vol*(Ainv(2,3)*Ainv(2,1) + Ainv(3,3)*Ainv(3,1)); % p3 mit p1 - selbe wie p1 mit p3
    B(T(i,3),T(i,2)) = B(T(i,3),T(i,2)) + tau*vol*(Ainv(2,3)*Ainv(2,2) + Ainv(3,3)*Ainv(3,2)); % p3 mit p2 - selbe wie p2 mit p3
    B(T(i,3),T(i,3)) = B(T(i,3),T(i,3)) + tau*vol*(Ainv(2,3)^2 + Ainv(3,3)^2); % p3 mit p3
    
    % 2. Teil der Bilinearform, d.h. die Integrale über phi_i*phi_j. Wir gehen wie in A 6.3 vor:
    for j = 1:3
        for k = 1:3
            % Baue zuerst das Produkt phi_i*phi_j als Funktion - jede
            % Spalte von Ainv beinhaltet die Polynomkoeffizienten einer
            % Hutfunktion.
            g = @(x,y) Ainv(1,k)*Ainv(1,j) + Ainv(2,k)*Ainv(2,j).*x.^2 + Ainv(3,k)*Ainv(3,j).*y.^2 + ...
                        Ainv(1,k)*Ainv(2,j).*x + Ainv(2,k)*Ainv(1,j).*x + ...
                        Ainv(1,k)*Ainv(3,j).*y + Ainv(3,k)*Ainv(1,j).*y + Ainv(2,k)*Ainv(3,j).*x.*y + Ainv(3,k)*Ainv(2,j).*x.*y;
            % werte nun mit eval_int das Integral aus: transofrmiere f auf
            % das Referenzdreieck und dort wird eine exakte Quadraturformel
            % verwendet. 2*vol entspricht der Trafodeterminante
            val = 2*vol*eval_int(g,P(T(i,:),:));
            % Anteil zur Steifigkeitsmatrix
            B(T(i,j),T(i,k)) = B(T(i,j),T(i,k)) + val; 
            % Entsprechenden Anteil aber auch speichern, da für den
            % Zeitabhängigen Teil von rechter Seite relevant, vgl. Theorie
            f_val((i-1)*3+j,k) = val; % könnte auch über (T(i,j),k) indizieren, ist aber nicht nötig.
        end;
    end;
    
    % 1. Teil der rechten Seite: zeitunabhängige Anteil der rechten Seite y
    % zuerst wertet man f in allen Punkten des Dreiecks aus: wenn alle
    % Punkte innerhlab des Teilgebietes sind (f > 0, bzw. f =1), so ist das Dreieck
    % drin und der Anteil des Dreiecks zur rechten Seite y wird bestimmt
    Punkte = P(T(i,:),:);
    Werte = f(Punkte);
    if isequal(Werte > 0, ones(3,1)) % Das Dreieck ist im Teilgebiet
        % der Flächeninhalt ist der Tetraeder-Inhalt, vgl. Theorieaufgabe
        y_part_one(T(i,1)) = y_part_one(T(i,1)) + tau/6*abs((P(T(i,2),1)-P(T(i,1),1))*(P(T(i,3),2)-P(T(i,1),2)) - (P(T(i,2),2)-P(T(i,1),2))*(P(T(i,3),1)-P(T(i,1),1)));
        y_part_one(T(i,2)) = y_part_one(T(i,2)) + tau/6*abs((P(T(i,1),1)-P(T(i,2),1))*(P(T(i,3),2)-P(T(i,2),2)) - (P(T(i,1),2)-P(T(i,2),2))*(P(T(i,3),1)-P(T(i,2),1)));
        y_part_one(T(i,3)) = y_part_one(T(i,3)) + tau/6*abs((P(T(i,1),1)-P(T(i,3),1))*(P(T(i,2),2)-P(T(i,3),2)) - (P(T(i,1),2)-P(T(i,3),2))*(P(T(i,2),1)-P(T(i,3),1)));
    end;
end;

% Randbehandlung
B(Dirichlet,:) = 0;
B(:,Dirichlet) = 0;
B = B + sparse(Dirichlet, Dirichlet, ones(size(Dirichlet)), n,n);
% einmalige LR-Zerlegung von B -> später dann BR\(BL\y) statt B\y in jedem
% Zeitschritt -> lohnt sich bie genügend Zeitschritten.
[BL,BR] = lu(B);

% Nun die Zeitevolution: in jedem Zeitschritt muss der 2. Teil der rechten
% Seite bestimmt werden und dann das LGS gelöst werden.
u = zeros(length(B),nt+1); % vgl. Anfangswerte in erster Spalte
for t = 1:nt
    y = y_part_one; % rechte Seite stets mit dem 1. part initialisieren    
    % baue den 2. Teil der rechten Seite - mit der Vorberechnung in f_val
    % ist es nun auf jedem Dreieck nur noch eine Matrix-Vektor-Multiplikation
    for i = 1:length(T)
        y(T(i,:)) = y(T(i,:)) + f_val((i-1)*3+1:(i-1)*3+3,:)*u(T(i,:),t);
    end;
    y(Dirichlet) = 0; % die 0-Randwerte müssen hier nochmals gesetzt werden, da Sie nicht aus den Koeffizienten von u(:,t) vererbt werden.
    % LGS lösen via LR-Zerlegung
    u(:,t+1) = BR\(BL\y);
end;

disp('Starte plot!')

% plot
g = figure;
for t = 1:nt+1
    % Es wird nur die Punktwolke geplottet (also die Werte an den Knoten).
    % Es fehlt die Information über die Dreiecke, ist aber schnell
   plot3(P(:,1),P(:,2),u(:,t),'.'); 
    axis([-4 4 -4 4 0 0.8]);    
    pause(0.01)
    
% Nachfolgend eine plot-variante in der jedes Dreieck in 3d nach und nach geplottet wird. Akkurat aber langsam.    
%     % nach und nach jedes Dreieck plotten. Dabei werden die Werte von u über die Farbkodierung wiedergegeben
    % hold on
    % for i = 1:length(T)
         % [1 2 3], da es immer Dreiecke sind und die Knoten sind entsprechen in P gespeichert.
    %     patch('Faces',[1 2 3],'Vertices',[P(T(i,:),:),u(T(i,:),t)],'FaceVertexCData',u(T(i,:),t),'EdgeColor','none','FaceColor','interp');
    % end;
    % view(3)
  %   pause
end;


function val = eval_int(f, points)
% - f ist die Funktion welche über dem Dreieck gegeben durch points integriert werden soll
% - points enthält die 3 Eckpunkte des Dreiecks
% - Anmerkung: die Trafo-Determinante beträgt 2*vol und wird außerhalb der Funktion dazumultipliziert.

% extrahiere die eckpunkte 
a = points(1,:);
b = points(2,:);
c = points(3,:);
% definiere die transformierte Funktion. Es ist f_transformiert =
% f(Phi(x)), wobei Phi das Referenzdreieck auf das Dreieck gegeben durch points abbildet.
f_transformiert = @(u,v) f(a(1)+u*(b(1)-a(1))+v*(c(1)-a(1)), a(2)+u*(b(2)-a(2))+v*(c(2)-a(2)));
% Damit kann man nun entweder via integral2 über das Refernzdreieck
% integrieren (auskommentiert), oder man verwendet eine in P_2 exakte
% Quadraturformel für das Referenzdreieck
%fy = @(x) 1-x;
%r = integral2(f_transformiert, 0,1, 0,fy);
val = 1/6*(f_transformiert(0.5,0.5) + f_transformiert(0.5,0) + f_transformiert(0,0.5)); % Quadraturformel


function I = Inhalt(p1,p2,p3)
% Flächeninhalt eines Dreiecks über die Determinante
a = p2-p1;
b = p3-p1;
I = 1/2*abs(a(1)*b(2) - a(2)*b(1));
