function eledata=element_output_xy(ielem,nodeCoordinates,...
    elementNodes,displacement,propMat,xi)
% xi is a vector containing the natural coordinates [xi eta] 
% where stress and strain will be computed

[numelem nelemnod]=size(elementNodes);

% Definition of hook law for plain stress
c1=propMat.E/(1-propMat.nu*propMat.nu);
c2=propMat.nu*c1;
G=propMat.E/(2*(1+propMat.nu));
C=[c1 c2 0;c2 c1 0; 0 0 G];

eledata.pos=zeros(1,2);
eledata.strain=zeros(1,3);
eledata.stress=zeros(1,3);

%Obtain strain at gauss point
xnod(1,1)=nodeCoordinates(elementNodes(ielem,1),1);
xnod(2,1)=nodeCoordinates(elementNodes(ielem,1),2);
xnod(3,1)=nodeCoordinates(elementNodes(ielem,2),1);
xnod(4,1)=nodeCoordinates(elementNodes(ielem,2),2);
xnod(5,1)=nodeCoordinates(elementNodes(ielem,3),1);
xnod(6,1)=nodeCoordinates(elementNodes(ielem,3),2);
xnod(7,1)=nodeCoordinates(elementNodes(ielem,4),1);
xnod(8,1)=nodeCoordinates(elementNodes(ielem,4),2);

unod(1,1)=displacement(elementNodes(ielem,1)*2-1);
unod(2,1)=displacement(elementNodes(ielem,1)*2);
unod(3,1)=displacement(elementNodes(ielem,2)*2-1);
unod(4,1)=displacement(elementNodes(ielem,2)*2);
unod(5,1)=displacement(elementNodes(ielem,3)*2-1);
unod(6,1)=displacement(elementNodes(ielem,3)*2);
unod(7,1)=displacement(elementNodes(ielem,4)*2-1);
unod(8,1)=displacement(elementNodes(ielem,4)*2);
B=Bmatrix_2D4N(xi,xnod);

% Computation of strain
epsilon=B*unod;

% Computation of stress
sigma=C*epsilon;

% Position in real coordinates (X,Y)
xpos=shapefunction_2D4N(xi,xnod);
eledata.pos=xpos';
eledata.strain=epsilon';
eledata.stress=sigma';
end