function A = assembDiffusion (eps, tri, P)

S1 =0.5* [ 1 -1 0; -1  1 0; 0 0 0];
S2 =0.5* [ 2 -1 -1; -1  0  1; -1  1  0];
S3 =0.5* [ 1  0 -1;  0  0  0; -1  0  1];

NuDoF = length(P);
A = sparse(NuDoF,NuDoF);
NuTri = length(tri);

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

B(:,1,:)=P(:,tri(:,2)) - P(:,tri(:,1));
B(:,2,:)=P(:,tri(:,3)) - P(:,tri(:,1));
Bdet = B(1,1,:).*B(2,2,:) - B(1,2,:).*B(2,1,:);
Bdet=abs(squeeze(Bdet));

gamma1 =   B(1,2,:).*B(1,2,:) + B(2,2,:).*B(2,2,:);
gamma2 = -(B(1,1,:).*B(1,2,:) + B(2,1,:).*B(2,2,:));
gamma3 =   B(2,1,:).*B(2,1,:) + B(1,1,:).*B(1,1,:);
gamma1=squeeze(gamma1)./Bdet;
gamma2=squeeze(gamma2)./Bdet;
gamma3=squeeze(gamma3)./Bdet;

for x=1:3
  for y=1:3
    temp = gamma1*S1(x,y)+gamma2*S2(x,y)+gamma3*S3(x,y);
    temp = temp.*eps;
    tempMatrix = sparse(tri(:,x), tri(:,y),temp,NuDoF,NuDoF);
    A = A + tempMatrix;
  end
end
