function  [deltaT ] = computeSDCoefficient(eps, k, mesh)

if length(eps) == 1
  eps = eps * ones(1,mesh.NuTri);
elseif length(eps) == mesh.NuDoF
  temp = (eps(mesh.T(:,1)) + eps(mesh.T(:,2)) + eps(mesh.T(:,3)))./3;
  eps = temp(:)';
  clear 'temp';
end

if  length(k) == mesh.NuDoF
  temp1 = (k(1,mesh.T(:,1)) + k(1,mesh.T(:,2)) + k(1,mesh.T(:,3)))./3;
  temp2 = (k(2,mesh.T(:,1)) + k(2,mesh.T(:,2)) + k(2,mesh.T(:,3)))./3;
  clear 'k';
  k(1,:) = temp1;
  k(2,:) = temp2;
end

firstNode  = mesh.edges(:,1);
secondNode = mesh.edges(:,2);
h = (mesh.P(1,firstNode)-mesh.P(1,secondNode) ).^2;
h = h + ( mesh.P(2,firstNode)-mesh.P(2,secondNode) ).^2;
h = sqrt(h);
hMax =  max(h(mesh.triangle2edge(1:mesh.NuTri,:))');

kNorm = sqrt(k(1,:).^2 + k(2,:).^2);

pecletZahl = 0.5*hMax.*kNorm./eps;
index =find(pecletZahl>1);
deltaT=zeros(size(pecletZahl));
deltaT(index) = hMax(index)./kNorm(index);
