%%Derivation of Linear triangular element matrices
%%This program does the algebra to get shape function (N)
%%and kinematic (B) matrices
%%Written by Sridhar Krishnaswamy 02/14/98
%
%
%%This program uses the Symbolic Toolbox of Matlab
%
%xy-coordinates defined as symbols
x=sym('x');
y=sym('y');
%Nodal Coordinates defined as symbols
x1=sym('X1'); y1=sym('Y1');
x2=sym('X2'); y2=sym('Y2');
x3=sym('X3'); y3=sym('Y3');
%
%"Rayleigh-Ritz" shape function is: Ux=a+bx+cy and similar for Uy
%Swap a,b,c in terms of nodal displacements U1,U2,U3 through
%{U1 U2 U3}'=Hmat*[a b c]' where: 
Hmat=[
1 x1 y1;
1 x2 y2;
1 x3 y3
]
%
% determinant of Hmat
disp('The determinant detH of the Hmatrix is:')
detH=det(Hmat)    
%Inverse of Hmat without the division by determinant to keep algebra not too messy
NNmat=detH.*inv(Hmat); 
%expose the following commented lines if you want to print each slot on a separate line;
%NN11=NNmat(1,1)
%NN12=NNmat(1,2)
%NN13=NNmat(1,3)
%NN21=NNmat(2,1)
%NN22=NNmat(2,2)
%NN23=NNmat(2,3)
%NN31=NNmat(3,1)
%NN32=NNmat(3,2)
%NN33=NNmat(3,3)
%
%The function Ux=(1/detH)[N1 N2 N3]{U1 U2 U3}' where:
N1=NNmat(1,1)+NNmat(2,1)*x+NNmat(3,1)*y;
N2=NNmat(1,2)+NNmat(2,2)*x+NNmat(3,2)*y;
N3=NNmat(1,3)+NNmat(2,3)*x+NNmat(3,3)*y;
%
%Using the same linear shape function for both Ux and Uy displacements:
%{Ux Uy}'=[Shape]{U1 V1 U2 V2 U3 V3}'
disp('The shape function matrix for the linear triangular element is:')
Shape=[
N1 0  N2 0  N3 0;
0  N1 0  N2 0  N3]
disp('The above matrix must be scalar multiplied by (1/detH)')
%Compute the B-matrix without (1/detH):
%Strain={epsxx epsyy gamxy}'={dUx/dx dUy/dy (dUy/dx+dUx/dy)}'
%Strain=[Bmat]{U1 V1 U2 V2 U3 V3}'
%Bmat=[	dN1/dx 	0 		dN2/dx	 0 		dN3/dx	0;
%      	0		dN1/dy	0		dN2/dy	0		dN3/dy;
%		dN1/dy	dN1/dx	dN2/dy	dN2/dx	dN3/dy	dN3/dx]
dN1dx=diff(N1,x); %this symbolically differentiates 'N1' wrt 'x' once
dN2dx=diff(N2,x);
dN3dx=diff(N3,x);
dN1dy=diff(N1,y);
dN2dy=diff(N2,y);
dN3dy=diff(N3,y);
%print Bmat
disp('The B-matrix for the linear triangular element is:')
Bmat=[	dN1dx 	0 		dN2dx	 0 		dN3dx	0;
      	0		dN1dy	0		dN2dy	0		dN3dy;
		dN1dy	dN1dx	dN2dy	dN2dx	dN3dy	dN3dx]
disp('The above matrix must be scalar multiplied by (1/detH)')

%endprogram