function Blatt_3

close all % schließe alle Grafikfenster
warning off % schalte Warnungen aus
addpath('~/Desktop/Lehre/NumPDGL_WS1617/FEM_Code_NPDGL/') % hier wird der Pfad angegeben, in welchem das 'distmesh'-Paket liegt
addpath('~/Desktop/Lehre/NumPDGL_WS1617/FEM_Code_NPDGL/') % hier wird der Pfad angegeben, in welchem das 'plotmesh'-Paket liegt

% Verschiedene Gebietsfunktionen
% Funktion für den Kreis um (5,4) mit Radius 1/sqrt(2)
f1 = @(p) sqrt((p(:,1) -5).^2 + (p(:,2) - 4).^2) - 1/sqrt(2);
% Funktion für ein Rechteck zwischen den Punkten p1 und p2
%   Anmerkung: In den 4 Flächen, welche jeweils von den 4 Ecken ausgehen
%   ist der eigentliche gerichtete Abstand die 2-norm zu jeder jeweiligen
%   Ecke. Dies ist hier nicht berücksichtig. Im Gegenzug geben wir immer
%   die Ecken als feste Punkte für das zu erstellende Gitter vor, vgl.
%   Dokumentation von distmesh2d.
p1 = [0,0]; p2 = [6,3];
f2 = @(p) -min(min(min(-p1(2) + p(:,2) , p2(2) - p(:,2)), -p1(1) + p(:,1)), p2(1) - p(:,1));
pfix1 = [p1;[0,3];[6,0];p2]; % für das zu f2 gehörende Rechteck
% Funktion für ein Rechteck zwischen den Punkten p3 und p4
p3 = [4,3]; p4 = [6,5];
f3 = @(p) -min(min(min(-p3(2) + p(:,2) , p4(2) - p(:,2)), -p3(1) + p(:,1)), p4(1) - p(:,1));
pfix2 = [p3;[6,3];[4,5];p4]; % für das zu f3 gehörende Rechteck

% Die Gebiete aus der Aufgabe
g1 = f1;
g2 = f2;
g3 = @(p) max(f3(p), -f1(p)); % f3 ohne f1
g4 = @(p) min(f2(p), f3(p)); % f2 vereinigt f3
pfix3 = [p1;[6,0];[0,3];p3;p4;[4,5]]; % Eckpunkte für die Vereinigung der beiden Rechtecke
%   Anmerkung: g4 gibt theoretisch nur dann den korrekten Abstand an, wenn
%   die beiden Gebiete, welche vereinigt werden disjunkt sind. Hier teilen
%   sich die beiden Rechtecke aber eine Kante, sodass dieser Fall
%   eigentlich noch berücksichtigt werden sollte. Dies füht in dieser
%   Variante dazu dass kein Punkt auf dieser Kante liegt.
% Die Gitter-Erstellung und der plot
% 1. Gebiet
figure;
[P1,T1]=distmesh2d(g1,@huniform,0.05,[2,2;6,6],[]);
plot_mesh(P1,T1);
disp('first domain done')
pause(2)
% 2. Gebiet
figure;
[P2,T2]=distmesh2d(g2,@huniform,0.1,[-1,-1;7,7],pfix1);
plot_mesh(P2,T2);
disp('second domain done')
pause(2)
% 3. Gebiet
figure;
[P3,T3]=distmesh2d(g3,@huniform,0.05,[2,2;6,6],pfix2);
plot_mesh(P3,T3);
disp('third domain done')
pause(2)
% 4. Gebiet
figure;
[P4,T4]=distmesh2d(g4,@huniform,0.1,[-1,-1;7,7],[pfix1;pfix3]);
plot_mesh(P4,T4);
disp('last domain done')
pause(2)

% Die Inhaltsberechnung

% Für die vier Gebiete aus der Aufgabe
I1 = Inhalt_Triangulierung(P1,T1);
I2 = Inhalt_Triangulierung(P2,T2);
I3 = Inhalt_Triangulierung(P3,T3);
I4 = Inhalt_Triangulierung(P4,T4);
% Für das Herz
load('grid.mat')
I5 = Inhalt_Triangulierung(p,t);
% Eine schöne Ausgabe
format long % mehr Nachkommastellen anzeigen lassen
fprintf('Inhalt des 1. Gebietes = %f. Exakter Inhalt = %f.\n' ,I1, pi/2)
fprintf('Inhalt des 2. Gebietes = %f. Exakter Inhalt = %f.\n' ,I2, 18)
fprintf('Inhalt des 3. Gebietes = %f. Exakter Inhalt = %f.\n' ,I3, 4-pi/2)
fprintf('Inhalt des 4. Gebietes = %f. Exakter Inhalt = %f.\n' ,I4, 18+4)
fprintf('Inhalt des Herzes = %f.\n' ,I5)

function I = Inhalt(p1,p2,p3)
% Über die Determinante
a = p2-p1;
b = p3-p1;
I = 1/2*abs(a(1)*b(2) - a(2)*b(1));

% Alternative über den Satz des Heron
% % Seitenlängen bestimmen
% a = norm(p2-p1);
% b = norm(p3-p1);
% c = norm(p3-p2);
% s = (a+b+c)/2; % halber Umfang des Dreiecks
% % Satz des Heron
% I = sqrt(s*(s-a)*(s-b)*(s-c));

function I = Inhalt_Triangulierung(P,T)

% Hilfsfunktion, welche mit Hilfe der Inhaltsfunktion den Inhalt der
% Triangulierung berechnet. Einfache for-Schleife über die Dreiecke.
I = 0;
for i = 1:length(T)
    I = I + Inhalt(P(T(i,1),:),P(T(i,2),:),P(T(i,3),:));
end
