function [divU,epsU] = diveps(x1,x2,x3,U1,U2,U3)

% [divU,epsU] = diveps(x1,x2,x3,U1,U2,U3)
%           calculates 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.

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