function Aloc = aloc_el(lambda,mu,xinit,e)

% Aloc = aloc_el(lambda,mu,xinit,e)
%             calculates a local contribution to the stiffness matrix
%             that corresponds to the bilinear form:
%              a(v,u) = lambda * \int_{K^e} div(v(x))*div(u(x)) dx +
%                     + mu * \int_{K^e} \sum_{k,l=1}^3 epsilon_{k,l}(v(x))*
%                                                      epsilon_{k,l}(u(x)) dx,
%             where eps_{k,l}(v) = 1/2 * ( dv_k/dx_l + dv_l/dx_k ).

x1 = xinit(:,e(1));
x2 = xinit(:,e(2));
x3 = xinit(:,e(3));

R = [x2-x1 x3-x1];
r = x1;

B = inv(R)'*[-1 1 0;-1 0 1];

Aloc = zeros(6,6);

for i1=1:3
    for j1=1:2
        for i2=1:3
            for j2=1:2
                index1 = (i1-1)*2 + j1;
                index2 = (i2-1)*2 + j2;
                
                Alambda = B(j1,i1)*lambda*B(j2,i2);
                
                Amu = 0;
                for k=1:2
                    for l=1:2
                        Bmu1 = 0;
                        Bmu2 = 0;
                        
                        if k==j1
                            Bmu1 = Bmu1 + 0.5*B(l,i1);
                        end
                        if l==j1
                            Bmu1 = Bmu1 + 0.5*B(k,i1);
                        end
                        if k==j2
                            Bmu2 = Bmu2 + 0.5*B(l,i2);
                        end
                        if l==j2
                            Bmu2 = Bmu2 + 0.5*B(k,i2);
                        end
                        
                        Amu = Amu + Bmu1*mu*Bmu2;
                    end
                end
                
                Aloc(index1,index2) = Aloc(index1,index2) + Alambda + Amu;
            end
        end
    end
end

Aloc = Aloc * (abs(det(R)) / 2);
