function [B,y,u]=Blatt5

clear all
close all

% Der folgende auskommentierte part fehlt in der online Version
% Der Code zur Erstellung der Baum-Triangulierung, sowie das finden der
% Dirichlet-Indizes.
% addpath('~/Desktop/Lehre/NumPDGL_WS1617/FEM_Code_NPDGL/') % hier wird der Pfad angegeben, in welchem das 'distmesh'-Paket liegt
% 
% % Polyeder für den Baum
% pfix = [[3,0];[3,2];[0,2];[3,5];[1,5];[3,7];[2,7];[4,9];[6,7];[5,7];[7,5];[5,5];[8,2];[5,2];[5,0];[3,0]];
% [P,T] = distmesh2d(@dpoly,@huniform,0.1,[0,0;9,9],pfix,pfix);
% 
% Dirichlet = find(abs(dpoly(P,pfix)) < 1e-3);

load('grid2.mat')

% Function handle, welcher das Teilgebiet K realisiert. Zu einem Punkt P
% ist f(P) <= 5*eps wenn P in K ist und > 5*eps wenn er außerhalb von K ist
f = @(p) min(sqrt((p(:,1) -4).^2 + (p(:,2) - 8).^2) - 0.5,eps) ...
        +min(sqrt((p(:,1) -3).^2 + (p(:,2) - 6).^2) - 0.5,eps) ...
        +min(sqrt((p(:,1) -4.5).^2 + (p(:,2) - 4).^2) - 0.5,eps) ...
        +min(sqrt((p(:,1) -2.5).^2 + (p(:,2) - 3).^2) - 0.5,eps)...
        +min(sqrt((p(:,1) -5).^2 + (p(:,2) - 6).^2) - 0.25,eps)...
        +min(sqrt((p(:,1) -6.5).^2 + (p(:,2) - 2.5).^2) - 0.25,eps);
% Initialisierung
n = length(P); % Anzahl Diskretisierungspunkte
B = spalloc(n,n,5*n); % Matrix B als sparse Matrix
y = zeros(n,1); % Rechte Seite -> kann später noch sparse gemacht werden, wobei die Anzahl der nichtnull-Einträge klar überwiegt

for i = 1:length(T)
    % Determinante aus dem LGS zur Bestimmung der Polynome p1,p2,p3
    % (jeweils 1 an einem der Eckpunkt und 0 in den anderen). Die Matrix
    % und damit die Determinante bleibt immer gleich, nur die rechte Seite
    % und damit die Einträge des Gradienten ändern sich, vgl. Theorieaufgabe
    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);
    % Flächeninahlt des Dreiecks
    vol = Inhalt(P(T(i,1),:),P(T(i,2),:),P(T(i,3),:));
    % die 6 Werte der Gradienten
    p1_x = 1/det*(P(T(i,2),2) - P(T(i,3),2));
    p1_y = 1/det*(P(T(i,3),1) - P(T(i,2),1));
    p2_x = 1/det*(P(T(i,3),2) - P(T(i,1),2));
    p2_y = 1/det*(P(T(i,1),1) - P(T(i,3),1));
    p3_x = 1/det*(P(T(i,1),2) - P(T(i,2),2));
    p3_y = 1/det*(P(T(i,2),1) - P(T(i,1),1));
    % Damit elementweises aufstellen der Matrix B
    B(T(i,1),T(i,1)) = B(T(i,1),T(i,1)) + vol*(p1_x^2 + p1_y^2); % p1 mit p1
    B(T(i,1),T(i,2)) = B(T(i,1),T(i,2)) + vol*(p1_x*p2_x + p1_y*p2_y); % p1 mit p2
    B(T(i,1),T(i,3)) = B(T(i,1),T(i,3)) + vol*(p1_x*p3_x + p1_y*p3_y); % p1 mit p3
    B(T(i,2),T(i,1)) = B(T(i,2),T(i,1)) + vol*(p2_x*p1_x + p2_y*p1_y); % p2 mit p1
    B(T(i,2),T(i,2)) = B(T(i,2),T(i,2)) + vol*(p2_x^2 + p2_y^2); % p2 mit p2
    B(T(i,2),T(i,3)) = B(T(i,2),T(i,3)) + vol*(p2_x*p3_x + p2_y*p3_y); % p2 mit p3
    B(T(i,3),T(i,1)) = B(T(i,3),T(i,1)) + vol*(p3_x*p1_x + p3_y*p1_y); % p3 mit p1
    B(T(i,3),T(i,2)) = B(T(i,3),T(i,2)) + vol*(p3_x*p2_x + p3_y*p2_y); % p3 mit p2
    B(T(i,3),T(i,3)) = B(T(i,3),T(i,3)) + vol*(p3_x^2 + p3_y^2); % p3 mit p3
    
    % nun noch die rechte Seite y
    % zuerst wertet man f in allen Punkten des Dreiecks aus: wenn alle
    % Punkte innerhlab des Teilgebietes sind (f <= 0), 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 <= 5*eps, ones(3,1)) % Das Dreieck ist im Teilgebiet
        % der Flächeninhalt ist der Tetraeder-Inhalt, vgl. Theorieaufgabe
        y(T(i,1)) = y(T(i,1)) + 1/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(T(i,2)) = y(T(i,2)) + 1/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(T(i,3)) = y(T(i,3)) + 1/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, vgl. auch Kommentar in der Übung
y(Dirichlet) = 0;
B(Dirichlet,:) = 0;
B(:,Dirichlet) = 0;
B = B + sparse(Dirichlet, Dirichlet, ones(size(Dirichlet)), n,n);
% LGS lösen
u = B\y;

% plot
g = figure;
% 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,:),:),'FaceVertexCData',u(T(i,:)),'EdgeColor','none','FaceColor','interp');
end
% Wir erstellen noch eine weihnachtliche Farbpalette
r=(0:(1/50):1)'; % rote Farbkomponente
g=(0.3:(0.6/50):0.9)'; % grüne Farbkomponente
b=(0.05:(0.3/50):0.35)'; % blaue Farbkomponente
map=[r,g,b]; % Colormap bauen
colormap(map) % Colormap verwenden
hold off


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));