P = [0 1 1 0 0.5; 0 0 1 1 0.5]; L = [3 4; 4 1; 1 2; 2 3];
T = [1 2 5; 1 4 5; 2 3 5; 3 4 5];
GL = zeros(4,1); GT = zeros(4,1);
[ mesh ] = buildMesh( T,P,[],GT,L,GL );
for i=1:2
  torefine = ones(mesh.NuTri,1);
  [ mesh] = refineMesh( mesh, torefine );
end
fehler=inf;

for i=3:7
  torefine = ones(mesh.NuTri,1);
  [ mesh] = refineMesh( mesh, torefine );

  x = mesh.P(1,:); y = mesh.P(2,:); 
  xT = ( x(mesh.T(:,1)) + x(mesh.T(:,2)) + x(1,mesh.T(:,3)) )/3;
  yT = ( y(mesh.T(:,1)) + y(mesh.T(:,2)) + y(1,mesh.T(:,3)) )/3;

  uLsg = 1-(x-0.5).^2 - (y-0.5).^2;

  eps = 10^-6; r = 1;
  k=[]; k(1,:) = xT; k(2,:) = yT.^2;
  f = -eps*(-4) -2*(x-0.5).*x - 2*(y-0.5).*y.^2 + r*uLsg;

  [ A, b ] = assembPDE( eps, k, r, f, mesh); 
  [ A, b ] = applyDirichletBoundaryConditions(A,b, mesh, 0,uLsg);
  u=A\b;

  alterFehler = fehler; 
  R = assembReaction (r, mesh.T, mesh.P);
  fehler = sqrt((u-uLsg')'*R*(u-uLsg'));
  fprintf('Fehler: %e  Faktor, %f \n',fehler,alterFehler/fehler);
end
trisurf(mesh.T,x,y,abs(u-uLsg'),'EdgeColor','none');
colormap(parula); colorbar;
xlabel('x'); ylabel('y');zlabel('z');
axis square; zlim([0 , max(abs(u-uLsg'))]);
view(-20, 25);
