function Blatt12

close all
dbstop if error

% initialisieren
h = 1/20;
g = @(x,y) -x-y;
f = @(x,y) 100 .* (abs(x - 0.5) < 0.05);
p = @(x,y) (x.^2 + y.^2 - 1).^3 - x.^2 .* y.^3;
% Berechne Approximation
u_h = level_set(f,g,p,h);

% Gitter übers Gebiet
[X,Y] = meshgrid(-2:h:2, -2:h:2);
X= X(:);
Y = Y(:);

%% verschiedene plots

% Variante 1
figure;
plot3(X,Y,u_h,'o')
xlabel('x')
ylabel('y')
zlabel('z')

% Variante 2 und 3 - es ist \Omega = (-2,2)^2 -> 4/h +1!
k = 4/h +1;
U = reshape(u_h,k,k);
X = reshape(X,k,k);
Y = reshape(Y,k,k);

% surf
figure;
surf(X,Y,U);
shading interp

% imagesc
U(p(X,Y) >= 0) = NaN;
figure;
imagesc(U);

% imagesc plottet die Matrixwerte. Diese korrespondieren aber eigentlich zu
% den FD-Punkten, welche im Gebiet von links unten nach rechts oben
% angeordnet. Das Bild müsste also noch entsprechend gedreht werden.




function u_h = level_set(f,g,p,h)

% Schrittweite konvertieren - die Breite des Gebietes ist 4
k = 4/h + 1;
% Gitter übers Gebiet
[X,Y] = meshgrid(-2:h:2, -2:h:2);
X= X(:);
Y = Y(:);

% Baue rechte Seite
f_h = f(X,Y);

% Baue Systemmatrix wie in Bem 2.5 / Bsp 2.6
T = spdiags(repmat([-1 2 -1],k,1),[-1 0 1],k,k);
A_h = 1/h^2*(kron(speye(k),T) + kron(T,speye(k)));

% Setzte nun Gleichungen, die zu Punkten x außerhalb von \Omega (gegeben
% durch p) liegen auf u(x) = g(x)
for i = 1:k^2
    if p(X(i),Y(i)) >= 0 % Punkt nicht im Gebiet
        A_h(i,:) = 0;
        A_h(i,i) = 1;
        f_h(i) = g(X(i),Y(i)); % erzwinge u(x) = g(x)
    end;
end;
% Berechne Approximation
u_h = A_h\f_h;
        