% ke_beam_sym: member stiffness matrix for grillage beam in global coordinates
syms x1 y1 x2 y2 real
syms h b th E nu positive 
dx = x2-x1;
dy = y2-y1;
L = (dx^2+dy^2)^(0.5);
cx = dx/L; cy = dy/L; 
ecc = th/2 + h/2;
G = E / (2*(1 + nu));
% Saint-Venant torsion constant
beta = 1/3 - 0.21*b/h*(1 - b^4/12/h^4); 
IX = beta*h*b^3; 
% Axial rigidity
EA_L = E*b*h/L;
% Rotational rigidities
GIX_L = G*IX/L;
IY = b*h^3/12;
EIY_L = E*IY/L;
EIY_L2 = EIY_L/L;
EIY_L3 = EIY_L2/L;
IZ = h*b^3/12;
EIZ_L = E*IZ/L;
EIZ_L2 = EIZ_L/L;
EIZ_L3 = EIZ_L2/L;
% 12-by-12 matrix in local coordinates
% u1,v1,w1,theta1,phi1,psi1,u2,v2,w2,theta2,phi2,psi2
sm = [EA_L 0 0 0 0 0 -EA_L 0 0 0 0 0;
      0 12*EIZ_L3 0 0 0 6*EIZ_L2 0 -12*EIZ_L3 0 0 0 6*EIZ_L2;
      0 0 12*EIY_L3 0 -6*EIY_L2 0 0 0 -12*EIY_L3 0 -6*EIY_L2 0;
      0 0 0 GIX_L 0 0 0 0 0 -GIX_L 0 0;
      0 0 -6*EIY_L2 0 4*EIY_L 0 0 0 6*EIY_L2 0 2*EIY_L 0;
      0 6*EIZ_L2 0 0 0 4*EIZ_L 0 -6*EIZ_L2 0 0 0 2*EIZ_L;
      -EA_L 0 0 0 0 0 EA_L 0 0 0 0 0;
      0 -12*EIZ_L3 0 0 0 -6*EIZ_L2 0 12*EIZ_L3 0 0 0 -6*EIZ_L2;
      0 0 -12*EIY_L3 0 6*EIY_L2 0 0 0 12*EIY_L3 0 6*EIY_L2 0;
      0 0 0 -GIX_L 0 0 0 0 0 GIX_L 0 0;
      0 0 -6*EIY_L2 0 2*EIY_L 0 0 0 6*EIY_L2 0 4*EIY_L 0;
      0 6*EIZ_L2 0 0 0 2*EIZ_L 0 -6*EIZ_L2 0 0 0 4*EIZ_L];
% Coupling terms
Tnum = eye(12);
T = sym(Tnum);
T(1,5) = ecc; 
T(2,4) = -ecc;
T(7,11) = ecc;
T(8,10) = -ecc;
smT = T'*(sm*T);
% Rotation matrix
R = [cx cy 0 0 0 0;
     -cy cx 0 0 0 0;
     0 0 1 0 0 0;
     0 0 0 cx cy 0;
     0 0 0 -cy cx 0;
     0 0 0 0 0 1];
RT = [R sym(zeros(6));
      sym(zeros(6)) R];
% Element matrix
ke = RT'*(smT*RT);

matlabFunction(ke,"File","ke_beam_from_sym");

% Partial derivatives
dkedx1 = diff(ke,x1);
matlabFunction(dkedx1,"File","dkedx1");
dkedy1 = diff(ke,y1);
matlabFunction(dkedy1,"File","dkedy1");
dkedx2 = diff(ke,x2);
matlabFunction(dkedx2,"File","dkedx2");
dkedy2 = diff(ke,y2);
matlabFunction(dkedy2,"File","dkedy2");
dkedh = diff(ke,h);
matlabFunction(dkedh,"File","dkedh");
dkedb = diff(ke,b);
matlabFunction(dkedb,"File","dkedb");
dkedth = diff(ke,th);
matlabFunction(dkedth,"File","dkedth");
