% ke_plate_sym: plate stiffness matrix in global coordinates
syms x1 y1 x2 y2 x3 y3 real
syms E nu th positive
%% Integration of DKT element
% Formulation following Batoz, Bathe, Ho 1980
x12 = x1-x2; x23 = x2-x3; x31 = x3-x1;
y12 = y1-y2; y23 = y2-y3; y31 = y3-y1;
l12 = (x12^2+y12^2)^0.5;
l23 = (x23^2+y23^2)^0.5;
l31 = (x31^2+y31^2)^0.5;
p6 = -6*x12/l12^2; p4 = -6*x23/l23^2; p5 = -6*x31/l31^2;
q6 = 3*x12*y12/l12^2; q4 = 3*x23*y23/l23^2; q5 = 3*x31*y31/l31^2;
t6 = -6*y12/l12^2; t4 = -6*y23/l23^2; t5 = -6*y31/l31^2;
r6 = 3*y12^2/l12^2; r4 = 3*y23^2/l23^2; r5 = 3*y31^2/l31^2;
% Integration points
IP = sym([0.5 0;0 0.5;0.5 0.5]);
% Loop
ke_DKT = sym(zeros(9));
for ip = 1:3
    xsi = IP(ip,1);
    eta = IP(ip,2);
    Hxdxsi = [p6*(1-2*xsi) + (p5-p6)*eta;
        q6*(1-2*xsi) - (q5+q6)*eta;
        -4 + 6*(xsi+eta) + r6*(1-2*xsi) - eta*(r5+r6);
        -p6*(1-2*xsi) + eta*(p4+p6);
        q6*(1-2*xsi) - eta*(q6-q4);
        -2 + 6*xsi + r6*(1-2*xsi) + eta*(r4-r6);
        -eta*(p5+p4);
        eta*(q4-q5);
        -eta*(r5-r4)];
    Hydxsi = [t6*(1-2*xsi) + (t5-t6)*eta;
        1 + r6*(1-2*xsi) - (r5+r6)*eta;
        -q6*(1-2*xsi) + eta*(q5+q6);
        -t6*(1-2*xsi) + eta*(t4+t6);
        -1 + r6*(1-2*xsi) + eta*(r4-r6);
        -q6*(1-2*xsi) - eta*(q4-q6);
        -eta*(t4+t5);
        eta*(r4-r5);
        -eta*(q4-q5)];
    Hxdeta = [-p5*(1-2*eta) - xsi*(p6-p5);
        q5*(1-2*eta) - xsi*(q5+q6);
        -4 + 6*(xsi+eta) + r5*(1-2*eta) - xsi*(r5+r6);
        xsi*(p4+p6);
        xsi*(q4-q6);
        -xsi*(r6-r4);
        p5*(1-2*eta) - xsi*(p4+p5);
        q5*(1-2*eta) + xsi*(q4-q5);
        -2 + 6*eta + r5*(1-2*eta) + xsi*(r4-r5)];
    Hydeta = [-t5*(1-2*eta) - xsi*(t6-t5);
        1 + r5*(1-2*eta) - xsi*(r5+r6);
        -q5*(1-2*eta) + xsi*(q5+q6);
        xsi*(t4+t6);
        xsi*(r4-r6);
        -xsi*(q4-q6);
        t5*(1-2*eta) - xsi*(t4+t5);
        -1 + r5*(1-2*eta) + xsi*(r4-r5);
        -q5*(1-2*eta) - xsi*(q4-q5)];
    % Area coefficient
    A = (x31*y12-x12*y31)/2;
    % Strain-displacement matrix
    Bmat = 1/2/A*[y31*Hxdxsi' + y12*Hxdeta';
        -x31*Hydxsi' - x12*Hydeta';
        -x31*Hxdxsi' - x12*Hxdeta' + y31*Hydxsi' + y12*Hydeta'];
    % Elasticity tensor
    Db = E*th^3/12/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
    % Integration
    intgrnd = Bmat'*(Db*Bmat);
    ke_DKT = ke_DKT + 1/3*A*intgrnd;
end
% Add fictitous stiffness
ke_DKT_full = sym(zeros(12));
indices = [1,2,3,5,6,7,9,10,11];
ke_DKT_full(indices,indices) = ke_DKT;
ke_DKT_full(4,4) = 1/1000*ke_DKT(1,1);
ke_DKT_full(8,8) = 1/1000*ke_DKT(4,4);
ke_DKT_full(12,12) = 1/1000*ke_DKT(7,7);

clear ke_DKT
ke_DKT = ke_DKT_full;
clear ke_DKT_full

%% Integration of CST element
Bmat = 1/2/A*[y23 0 y31 0 y12 0;
    0 -x23 0 -x31 0 -x12;
    -x23 y23 -x31 y31 -x12 y12];
Ds = E/(1-nu^2)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2];
ke_CST = 1*A*th*Bmat'*(Ds*Bmat);

%% Numerical function
matlabFunction(ke_DKT,"File","ke_DKT_from_sym");
matlabFunction(ke_CST,"File","ke_CST_from_sym");

%% Partial derivatives
% DKT
dkeDKTdx1 = diff(ke_DKT,x1);
matlabFunction(dkeDKTdx1,"File","dkeDKTdx1");
dkeDKTdy1 = diff(ke_DKT,y1);
matlabFunction(dkeDKTdy1,"File","dkeDKTdy1");
dkeDKTdx2 = diff(ke_DKT,x2);
matlabFunction(dkeDKTdx2,"File","dkeDKTdx2");
dkeDKTdy2 = diff(ke_DKT,y2);
matlabFunction(dkeDKTdy2,"File","dkeDKTdy2");
dkeDKTdx3 = diff(ke_DKT,x3);
matlabFunction(dkeDKTdx3,"File","dkeDKTdx3");
dkeDKTdy3 = diff(ke_DKT,y3);
matlabFunction(dkeDKTdy3,"File","dkeDKTdy3");
dkeDKTdth = diff(ke_DKT,th);
matlabFunction(dkeDKTdth,"File","dkeDKTdth");
% CST
dkeCSTdx1 = diff(ke_CST,x1);
matlabFunction(dkeCSTdx1,"File","dkeCSTdx1");
dkeCSTdy1 = diff(ke_CST,y1);
matlabFunction(dkeCSTdy1,"File","dkeCSTdy1");
dkeCSTdx2 = diff(ke_CST,x2);
matlabFunction(dkeCSTdx2,"File","dkeCSTdx2");
dkeCSTdy2 = diff(ke_CST,y2);
matlabFunction(dkeCSTdy2,"File","dkeCSTdy2");
dkeCSTdx3 = diff(ke_CST,x3);
matlabFunction(dkeCSTdx3,"File","dkeCSTdx3");
dkeCSTdy3 = diff(ke_CST,y3);
matlabFunction(dkeCSTdy3,"File","dkeCSTdy3");
dkeCSTdth = diff(ke_CST,th);
matlabFunction(dkeCSTdth,"File","dkeCSTdth");