function A = assembReaction (r, tri, P)

r = r(:);
S4= [2, 1, 1; 1, 2, 1; 1, 1, 2]/24;

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(r) == 1
  temp=r;
  r = ones(NuTri,1);
  r = temp*r;
  clear 'temp';
elseif length(r) == NuDoF
  temp = (r(tri(:,1)) + r(tri(:,2)) + r(tri(:,3)))./3;
  r = temp;
  clear 'temp';
end

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));

for x=1:3
  for y=1:3
    temp = Bdet*S4(x,y);
    temp = temp.*r;
    tempMatrix = sparse(tri(:,x), tri(:,y),temp,NuDoF,NuDoF);
    A = A + tempMatrix;
  end
end