%function [k]=PlateQ4_2Dk(X, Y, t, E, Nu)
%clear all;                                   
X=[50 100 150 70];
Y=[60 70 180 120];
t =10;
E=20000;
Nu=0.3;

D=E/(1-Nu^2)*[1 Nu 0 0 0;Nu 1 0 0 0; 0 0 (1-Nu)/2 0 0; 0 0 0 (1-Nu)/2 0; 0 0 0 0 (1-Nu)/2];

syms r s z real;



%XiC=[-1 1 1 -1];
%EtaC=[-1 -1 1 1];

%Ni=[ ( (XiC(2)-Xi)*(EtaC(4)-Eta) )/( (XiC(2)-XiC(1))*(EtaC(4)-EtaC(1))),    ((XiC(1)-Xi)*(EtaC(3)-Eta))/( (XiC(1)-XiC(2))*(EtaC(3)-EtaC(2))),      ((XiC(4)-Xi)*(EtaC(2)-Eta))/( (XiC(4)-XiC(3))*(EtaC(2)-EtaC(3))),      ((XiC(3)-Xi)*(EtaC(1)-Eta))/( (XiC(3)-XiC(4))*(EtaC(1)-EtaC(4))) ];
Ni=[0.25*(1-r)*(1-s) 0.25*(1+r)*(1-s) 0.25*(1+r)*(1+s) 0.25*(1-r)*(1+s) ];

dNibydr=diff(Ni,r);
dNibyds=diff(Ni,s);


x=Ni*X';
y=Ni*Y';

J=[ diff(x,r), diff(y,r); diff(x,s), diff(y,s)];
Jinv=inv(J);

dNibydx=Jinv(1,1)*dNibydr+Jinv(1,2)*dNibyds;
dNibydy=Jinv(2,1)*dNibydr+Jinv(2,2)*dNibyds;


B1=[0     -z*dNibydx(1)         0           0      -z*dNibydx(2)        0             0      -z*dNibydx(3)       0           0      -z*dNibydx(4)         0 ; ...
   0           0        -z*dNibydy(1)      0            0          -z*dNibydy(2)     0           0         -z*dNibydy(3)    0            0         -z*dNibydy(4);   ...
   0      -z*dNibydy(1)   -z*dNibydx(1)    0      -z*dNibydy(2)    -z*dNibydx(2)     0      -z*dNibydy(3)  -z*dNibydx(3)    0      -z*dNibydy(4)   -z*dNibydx(4); ...
   dNibydy(1)  0            -Ni(1)     dNibydy(2)      0               -Ni(2)     dNibydy(3)     0            -Ni(3)     dNibydy(4)      0            -Ni(4); ...
   dNibydx(1)  -Ni(1)         0        dNibydx(2)      -Ni(2)          0          dNibydx(3)     -Ni(3)         0        dNibydx(4)      -Ni(4)         0];    

BDB=int(B'*D*B,z,-5,+5);

BTDBJ1=subs(BDB*det(J),{r, s},{ -1/sqrt(3), -1/sqrt(3)});
BTDBJ2=subs(BDB*det(J),{r, s},{ 1/sqrt(3), -1/sqrt(3)});
BTDBJ3=subs(BDB*det(J),{r, s},{ 1/sqrt(3), 1/sqrt(3)});
BTDBJ4=subs(BDB*det(J),{r, s},{ -1/sqrt(3), 1/sqrt(3)});
k=BTDBJ1+BTDBJ2+BTDBJ3+BTDBJ4;





   