function [l2,h1,linfty] = Fehler(P, T, u, ux, uy, uh)

% Berechnet den Fehler zur exakten Loesung
% Aufruf: [l2,h1] = Fehler(P, T, u, ux, uy, uh)
%
%   P:  Gitter-Punkte
%   T:  Dreiecke
%   u:  Matlab-Datei (Name als Zeichenkette)
%   ux: Matlab-Datei (Name als Zeichenkette)
%   uy: Matlab-Datei (Name als Zeichenkette)
%   uh: Diskrete Loesung

format long
% Anzahl der Elemente und der Eckpunkte
NEL = size(T,2);
NVT = size(P,2);

% initialisiere Variablen
Area = 0;
l2 = 0;
h1 = 0;
linfty = 0;

% calculate gradient of uh using PDEGRAD
[uhx, uhy] = pdegrad(P, T, uh);

% fuer alle Dreiecke
for i = 1:NEL,
% set lokale Variablen auf Null 
% locl2 = local L^2
% loch1 = local H^1
% locinf = local L^infty

  locl2 = 0;
  loch1 = 0;
  locinf = 0;

% berechne Eckpunkte des Dreiecks
  V1 = T(1,i);
  P1 = P(1:2,V1);
  x1 = P1(1);
  y1 = P1(2);

  V2 = T(2,i);
  P2 = P(1:2,V2);
  x2 = P2(1);
  y2 = P2(2);

  V3 = T(3,i);
  P3 = P(1:2,V3);
  x3 = P3(1);
  y3 = P3(2);

% berechne Inhalt des Dreiecks
  A = 0.5*((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1));
  Area = Area + A;

% Finite Element Approximation in Eckpunkten
  val1 =  uh(V1);
  val2 =  uh(V2);
  val3 =  uh(V3);
% Loesung in den Eckpunkten
  x = [x1 x2 x3];
  y = [y1 y2 y3];
  ue = feval(u,x,y);
% L^infty Fehler
  locinf = max(abs(ue(1)-val1),abs(ue(2)-val2));
  locinf = max(locinf,abs(ue(3)-val3));

% berechne die Werte der Loesung der der FE-Approximation
% in den Mittelpunkten der Kanten

% Mittelpunkte der Kanten
  e1x = 0.5*(x1+x2);
  e1y = 0.5*(y1+y2);
  e2x = 0.5*(x2+x3);
  e2y = 0.5*(y2+y3);
  e3x = 0.5*(x1+x3);
  e3y = 0.5*(y1+y3);

% Finite Element Approximation
  val1 = 0.5*(uh(V1)+uh(V2));
  val2 = 0.5*(uh(V2)+uh(V3));
  val3 = 0.5*(uh(V1)+uh(V3));

% Ableitungen der FE-Approximation (konstant)
  valx = uhx(i);
  valy = uhy(i);

  x = [e1x e2x e3x];
  y = [e1y e2y e3y];

% (exakte) Loesung
  ue = feval(u,x,y);
  uex = feval(ux,x,y);
  uey = feval(uy,x,y);

% berechne lokale Fehler und akkumuliere globale Fehler
  locl2 = A * sum( (ue-[val1 val2 val3]).^2 ) /3;
  l2 = l2 + locl2;

  loch1 = A * ( sum((uex-valx).^2) + sum((uey-valy).^2) ) /3;
  h1 = h1 + loch1;

% L^infty Fehler
%  locinf = max(locinf,abs(ue(1)-val1));
%  locinf = max(locinf,abs(ue(2)-val2));
%  locinf = max(locinf,abs(ue(3)-val3));

  if (locinf > linfty)
    linfty = locinf;
  end
end;

% ziehe Quadratwurzel
l2 = sqrt(l2);
h1 = sqrt(h1);

format loose
