function [graddivU,gradepsU] = graddiveps(x1,x2,x3,U1,U2,U3)

% [graddivU,gradepsU] = graddiveps(x1,x2,x3,U1,U2,U3)
%           calculates gradient of both the divergence and screw of the linear
%           function U given by its 2d values U1,U2,U3
%           at the corners x1,x2,x3 of a triangular element
%           with respect to U1, U2, and U3.

R = [x2-x1 , x3-x1];
invR = inv(R);

%dUdx1 = (U2-U1)*invR(1,1) + (U3-U1)*invR(2,1);
%dUdx2 = (U2-U1)*invR(1,2) + (U3-U1)*invR(2,2);

%divU = dUdx1(1) + dUdx2(2);
%curlU = 0.5 * (dUdx1(2) + dUdx2(1));
%epsU = [dUdx1(1) , curlU ; curlU , dUdx2(2) ];

graddivU = [-invR(1,1)-invR(2,1);-invR(1,2)-invR(2,2);invR(1,1);invR(1,2);invR(2,1);invR(2,2)];

gradepsU = zeros(2,12);
% depsU/dU11
gradepsU(:,1:2) = [-invR(1,1)-invR(2,1) 0.5*(-invR(1,2)-invR(2,2));0.5*(-invR(1,2)-invR(2,2)) 0];
% depsU/dU12
gradepsU(:,3:4) = [0 0.5*(-invR(1,1)-invR(2,1));0.5*(-invR(1,1)-invR(2,1)) -invR(1,2)-invR(2,2)];
% depsU/dU21
gradepsU(:,5:6) = [invR(1,1) 0.5*invR(1,2);0.5*invR(1,2) 0];
% depsU/dU22
gradepsU(:,7:8) = [0 0.5*invR(1,1);0.5*invR(1,1) invR(1,2)];
% depsU/dU31
gradepsU(:,9:10) = [invR(2,1) 0.5*invR(2,2);0.5*invR(2,2) 0];
% depsU/dU32
gradepsU(:,11:12) = [0 0.5*invR(2,1);0.5*invR(2,1) invR(2,2)];
