function w = QuadTriangleElementSresses(E,NU,t,x1,y1,x2,y2,x3,y3,p,u)
%QuadTriangleElementStresses   This function returns the element 
%                              stresses for a quadratic   
%                              triangular element with modulus of 
%                              elasticity E, Poisson's ratio NU,  
%                              thickness t, coordinates of the 
%                              node 1 (x1,y1), coordinates of 
%                              node 2 (x2,y2),  
%                              coordinates of node 3 
%                              (x3,y3), and element displacement
%                              vector u. Use p = 1 for cases of 
%                              plane stress, and p = 2 for cases 
%                              of plane strain.
%                              The size of the element stiffness 
%                              matrix is 12 x 12.
syms x y;
x4 = (x1 + x2)/2;
y4 = (y1 + y2)/2;
x5 = (x2 + x3)/2;
y5 = (y2 + y3)/2;
x6 = (x1 + x3)/2;
y6 = (y1 + y3)/2;
x21 = x2 - x1;
y21 = y2 - y1;
x23 = x2 - x3;
y23 = y2 - y3;
x46 = x4 - x6;
y46 = y4 - y6;
x13 = x1 - x3;
y13 = y1 - y3;
x16 = x1 - x6;
y16 = y1 - y6;
x31 = x3 - x1;
y31 = y3 - y1;
x54 = x5 - x4;
y54 = y5 - y4;
x24 = x2 - x4;
y24 = y2 - y4;
x56 = x5 - x6;
y56 = y5 - y6;
x36 = x3 - x6;
y36 = y3 - y6;
x41 = x4 - x1;
y41 = y4 - y1;
x43 = x4 - x3;
y43 = y4 - y3;
x51 = x5 - x1;
y51 = y5 - y1;
x61 = x6 - x1;
y61 = y6 - y1;
x63 = x6 - x3;
y63 = y6 - y3;
N1 = (x23*(y-y3)-y23*(x-x3))*(x46*(y-y6)-y46*(x-x6))/((x23*y13-y23*x13)*(x46*y16-y46*x16));
N2 = (x31*(y-y1)-y31*(x-x1))*(x54*(y-y4)-y54*(x-x4))/((x31*y21-y31*x21)*(x54*y24-y54*x24));
N3 = (x21*(y-y1)-y21*(x-x1))*(x56*(y-y6)-y56*(x-x6))/((x21*y31-y21*x31)*(x56*y36-y56*x36));
N4 = (x31*(y-y1)-y31*(x-x1))*(x23*(y-y3)-y23*(x-x3))/((x31*y41-y31*x41)*(x23*y43-y23*x43));
N5 = (x31*(y-y1)-y31*(x-x1))*(x21*(y-y1)-y21*(x-x1))/((x31*y51-y31*x51)*(x21*y51-y21*x51));
N6 = (x21*(y-y1)-y21*(x-x1))*(x23*(y-y3)-y23*(x-x3))/((x21*y61-y21*x61)*(x23*y63-y23*x63));
N1x = diff(N1,x);
N1y = diff(N1,y);
N2x = diff(N2,x);
N2y = diff(N2,y);
N3x = diff(N3,x);
N3y = diff(N3,y);
N4x = diff(N4,x);
N4y = diff(N4,y);
N5x = diff(N5,x);
N5y = diff(N5,y);
N6x = diff(N6,x);
N6y = diff(N6,y);
B = [N1x, 0, N2x, 0, N3x, 0, N4x, 0, N5x, 0, N6x, 0 ;
   0, N1y, 0, N2y, 0, N3y, 0, N4y, 0, N5y, 0, N6y;
   N1y, N1x, N2y, N2x, N3y, N3x, N4y, N4x, N5y, N5x, N6y, N6x];
if p == 1 
   D = (E/(1-NU*NU))*[1, NU, 0 ; NU, 1, 0 ; 0, 0, (1-NU)/2];
elseif p == 2
   D = (E/(1+NU)/(1-2*NU))*[1-NU, NU, 0 ; NU, 1-NU, 0 ; 0, 0, (1-2*NU)/2];
end
w = D*B*u
%
% We also calculate the stresses at the centroid of the element
%
xcent = (x1 + x2 + x3)/3;
ycent = (y1 + y2 + y3)/3;
wcent = subs(w, {x,y}, {xcent,ycent});
w = double(wcent);