function gradAloc_mu = gradAloc_mu(mu,x,e)

% gradAloc_mu = gradAloc_mu(mu,x,e)
%             It calculates gradient of a local contribution to the system matrix
%             that corresponds to the bilinear form:
%              a(v,u) = (1/mu) * \int_{K^e} grad(v(x))*grad(u(x)) dx,
%             with respect to nodal coordinates.

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

R = [x2-x1 x3-x1];
invR = inv(R);
detR = det(R);
graddetR = [x2(2)-x3(2) x3(2)-x1(2) x1(2)-x2(2);x3(1)-x2(1) x1(1)-x3(1) x2(1)-x1(1)]';
Rtilde = [x3(2)-x1(2) x1(2)-x2(2);x1(1)-x3(1) x2(1)-x1(1)];

D = abs(detR) * (1/mu);
B = invR' * [-1 1 0;-1 0 1];

gradD = (sign(detR) * (1/mu)) * graddetR;
    
dinvRdx11 = (1/power(detR,2)) * ( detR*[0 0;1 -1]' - graddetR(1,1)*Rtilde' );
dinvRdx12 = (1/power(detR,2)) * ( detR*[-1 1;0 0]' - graddetR(1,2)*Rtilde' );
dinvRdx21 = (1/power(detR,2)) * ( detR*[0 0;0 1]' - graddetR(2,1)*Rtilde' );
dinvRdx22 = (1/power(detR,2)) * ( detR*[0 -1;0 0]' - graddetR(2,2)*Rtilde' );
dinvRdx31 = (1/power(detR,2)) * ( detR*[0 0;-1 0]' - graddetR(3,1)*Rtilde' );
dinvRdx32 = (1/power(detR,2)) * ( detR*[1 0;0 0]' - graddetR(3,2)*Rtilde' );

dBdx11 = dinvRdx11'*[-1 1 0;-1 0 1];
dBdx12 = dinvRdx12'*[-1 1 0;-1 0 1];
dBdx21 = dinvRdx21'*[-1 1 0;-1 0 1];
dBdx22 = dinvRdx22'*[-1 1 0;-1 0 1];
dBdx31 = dinvRdx31'*[-1 1 0;-1 0 1];
dBdx32 = dinvRdx32'*[-1 1 0;-1 0 1];

gradAloc_mu = zeros(3,18);

gradAloc_mu(:,1:3) = 0.5 * ( dBdx11'*D*B + B'*gradD(1,1)*B + B'*D*dBdx11 );
gradAloc_mu(:,4:6) = 0.5 * ( dBdx12'*D*B + B'*gradD(1,2)*B + B'*D*dBdx12 );
gradAloc_mu(:,7:9) = 0.5 * ( dBdx21'*D*B + B'*gradD(2,1)*B + B'*D*dBdx21 );
gradAloc_mu(:,10:12) = 0.5 * ( dBdx22'*D*B + B'*gradD(2,2)*B + B'*D*dBdx22 );
gradAloc_mu(:,13:15) = 0.5 * ( dBdx31'*D*B + B'*gradD(3,1)*B + B'*D*dBdx31 );
gradAloc_mu(:,16:18) = 0.5 * ( dBdx32'*D*B + B'*gradD(3,2)*B + B'*D*dBdx32 );

%Aloc_mu = 0.5 * B' * D * B;
