function gradbloc_el = gradbloc_el(lambda,mu,graddivU,gradepsU,xinit,e)

% gradbloc_el = gradbloc_el(lambda,nu,graddivU,gradepsU,xinit,e)
%             calculates gradient w.r. to U1, U2, and U3
%             of a local element contribution to the
%             right-hand side vector that corresponds
%             to the following linear functional:
%              f(v) = - lambda * \int_{K^e} divU*div(v(x)) dx +
%                     - mu * \int_{K^e} \sum_{k,l=1}^3 epsU_{k,l}*eps_{k,l}(v(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];

gradbloc = zeros(6,6);

for i1=1:3
    for j1=1:2
        index1 = (i1-1)*2 + j1;
        
%        blambda = divU*lambda*B(j1,i1)*graddivU;
        gradblambda = lambda*B(j1,i1) * graddivU;

%        bmu = 0;
        gradbmu = zeros(6,1);
        for k=1:2
            for l=1:2
                Bmu1 = 0;

                if k==j1 
                    Bmu1 = Bmu1 + 0.5*B(l,i1);
                end
                if l==j1
                    Bmu1 = Bmu1 + 0.5*B(k,i1);
                end
                
%                bmu = bmu + epsU(k,l)*mu*Bmu1;
                gradbmu = gradbmu + mu*Bmu1 * [gradepsU(k,l);gradepsU(k,2+l);gradepsU(k,4+l);gradepsU(k,6+l);gradepsU(k,8+l);gradepsU(k,10+l)];
            end
        end

%        bloc(index1) = bloc(index1) + blambda + bmu;
        gradbloc(index1,:) = gradbloc(index1,:) + gradblambda' + gradbmu';

    end
end

%bloc_el = - bloc * (abs(det(R)) / 2);
gradbloc_el = - gradbloc * (abs(det(R)) / 2);
