clear;clc

nInterp = 20; % Interpolationsaufloesung
nPlot = 100;  % Plotaufloesung fuer v in L^2(Omega)

% Definiere Funktion mit Sprung.
v_links = @(x,y) x.^4-sqrt(y)-sin(x.*y.*2*pi);
v_rechts = @(x,y) -x.^4+sqrt(y)+sin(x.*y.*2*pi);
v = @(x,y) (x > 0.5).*v_rechts(x,y) + (x <= 0.5).*v_links(x,y);

% Erzeuge Triangulierung (FEM-Gitter).
xlin = linspace(0,1,nInterp+1);
[xx,yy] = meshgrid(xlin,xlin);
tri = delaunay(xx(:),yy(:));
x = [xx(:),yy(:)];

% Quadraturformel: exakt bis Polynomgrad 2
quadratur_P2.knoten = [
	4,1;
	1,4;
	1,1]/6;
quadratur_P2.gewichte = [1;1;1]/6;

% Quadraturformel: exakt bis Polynomgrad 3
quadratur_P3.knoten = [
    1/3,1/3;
    1/5,3/5;
    1/5,1/5;
    3/5,1/5];
quadratur_P3.gewichte = [-27/48;25/48;25/48;25/48]/2;

% Reserviere Speicher fuer Interpolation p in P1(Omega).
p = zeros(size(x,1),1);

% Bestimme Interpolationswerte (Pj v)(xj) =: pj(xj).
for j = 1:size(x,1)
    % Bestimme Dreiecke, die den Knoten xj enthalten.
    istNachbarDreieck = any(tri == j,2);
    nachbarDreiecke = find(istNachbarDreieck);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% Hier lokale Projektion (Pj v) = pj und p(xj) := pj(xj) berechnen.
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

% Erzeuge Gitter zum Plotten der Funktion v.
xlin = linspace(0,1,nPlot+1);
[xx,yy] = meshgrid(xlin,xlin);
tri_v = delaunay(xx(:),yy(:));
x_v = [xx(:),yy(:)];

% Werte die Funktion v im linken/rechten Teil aus.
veval_l = v_links( x_v(:,1),x_v(:,2));
veval_r = v_rechts(x_v(:,1),x_v(:,2));

% Bestimme alle Dreiecke, die im linken Bereich liegen.
x1 = x_v(:,1);
tri_l = all(x1(tri_v) <= 0.5,2);

figure

% Plotte Funktion v in linken Breich.
trisurf(tri_v(tri_l,:),x_v(:,1),x_v(:,2),veval_l,'FaceAlpha',0.7,'EdgeColor','none')
hold on

% Plotte Funktion v in rechten Breich.
trisurf(tri_v(not(tri_l),:),x_v(:,1),x_v(:,2),veval_r,'FaceAlpha',0.7,'EdgeColor','none')

% Plotte Interpolation p.
trisurf(tri,x(:,1),x(:,2),p,'FaceColor','none','EdgeColor','r','LineWidth',2)

% Plotfarben korrigieren.
caxis([min(min(veval_l(:)),min(veval_r(:))), ...
       max(max(veval_l(:)),max(veval_r(:)))])
