function [k]=Q8_3Dk(X, Y, Z, E, Nu)

%X=[5 12 13 7 5 12 13 7];
%Y=[6 8 11 12 6 8 11 12];
%Z=[-10 -10 -10 -10 10 10 10 10];
%E=2e11;
%Nu=0.25;
%Q8_3Dk(X,Y,Z,E,Nu)
%


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


syms Xi Eta Zeta real;



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

Ni=[ ( (XiC(2)-Xi)*(EtaC(4)-Eta)*(ZetaC(5)-Zeta ) )/( (XiC(2)-XiC(1))*(EtaC(4)-EtaC(1))*(ZetaC(5)-ZetaC(1))), ...
        ( (XiC(1)-Xi)*(EtaC(3)-Eta)*(ZetaC(6)-Zeta ) )/( (XiC(1)-XiC(2))*(EtaC(3)-EtaC(2))*(ZetaC(6)-ZetaC(2))), ...
        ( (XiC(4)-Xi)*(EtaC(2)-Eta)*(ZetaC(7)-Zeta ) )/( (XiC(4)-XiC(3))*(EtaC(2)-EtaC(3))*(ZetaC(7)-ZetaC(3))), ...
        ( (XiC(3)-Xi)*(EtaC(1)-Eta)*(ZetaC(8)-Zeta ) )/( (XiC(3)-XiC(4))*(EtaC(1)-EtaC(4))*(ZetaC(8)-ZetaC(4))), ...
        ( (XiC(6)-Xi)*(EtaC(8)-Eta)*(ZetaC(1)-Zeta ) )/( (XiC(6)-XiC(5))*(EtaC(8)-EtaC(5))*(ZetaC(1)-ZetaC(5))), ...
        ( (XiC(5)-Xi)*(EtaC(7)-Eta)*(ZetaC(2)-Zeta ) )/( (XiC(5)-XiC(6))*(EtaC(7)-EtaC(6))*(ZetaC(2)-ZetaC(6))), ...
        ( (XiC(8)-Xi)*(EtaC(6)-Eta)*(ZetaC(3)-Zeta ) )/( (XiC(8)-XiC(7))*(EtaC(6)-EtaC(7))*(ZetaC(3)-ZetaC(7))), ...
        ( (XiC(7)-Xi)*(EtaC(5)-Eta)*(ZetaC(4)-Zeta ) )/( (XiC(7)-XiC(8))*(EtaC(5)-EtaC(8))*(ZetaC(4)-ZetaC(8)))];
    
dNibydXi=diff(Ni,Xi);
dNibydEta=diff(Ni,Eta);
dNibydZeta=diff(Ni,Zeta);



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

J=[ diff(x,Xi), diff(y,Xi), diff(z,Xi) ; diff(x,Eta), diff(y,Eta), diff(z,Eta); diff(x,Zeta), diff(y,Zeta), diff(z,Zeta)];
Jinv=inv(J);

dNibydx=simplify(Jinv(1,1)*dNibydXi+Jinv(1,2)*dNibydEta+Jinv(1,3)*dNibydZeta);
dNibydy=simplify(Jinv(2,1)*dNibydXi+Jinv(2,2)*dNibydEta+Jinv(2,3)*dNibydZeta);
dNibydz=simplify(Jinv(3,1)*dNibydXi+Jinv(3,2)*dNibydEta+Jinv(3,3)*dNibydZeta);

B=[ dNibydx(1)         0                       0            dNibydx(2)           0                     0              dNibydx(3)             0                        0             dNibydx(4)          0                        0   ...
       dNibydx(5)          0                      0             dNibydx(6)           0                     0              dNibydx(7)             0                       0             dNibydx(8)           0                       0;   ...
             0               dNibydy(1)           0                     0             dNibydy(2)           0                       0               dNibydy(3)            0                     0               dNibydy(4)           0   ...
             0               dNibydy(5)           0                     0             dNibydy(6)           0                       0               dNibydy(7)            0                     0               dNibydy(8)           0;   ...      
             0                    0                dNibydz(1)           0                      0             dNibydz(2)            0                      0               dNibydz(3)           0                      0               dNibydz(4)  ...
             0                    0                dNibydz(5)           0                      0             dNibydz(6)            0                      0               dNibydz(7)           0                      0               dNibydz(8);   ...      
             0               dNibydz(1)   dNibydy(1)           0             dNibydz(2)   dNibydy(2)            0               dNibydz(3)   dNibydy(3)            0              dNibydz(4)    dNibydy(4) ...
             0               dNibydz(5)   dNibydy(5)           0             dNibydz(6)   dNibydy(2)            0               dNibydz(7)   dNibydy(2)            0              dNibydz(8)    dNibydy(2);   ...                   
       dNibydx(1)        0                dNibydz(1)   dNibydx(2)           0             dNibydz(2)    dNibydx(3)           0               dNibydz(3)   dNibydx(4)            0               dNibydz(4)  ...
       dNibydx(5)        0                dNibydz(5)   dNibydx(6)           0             dNibydz(6)    dNibydx(7)           0               dNibydz(7)   dNibydx(8)            0               dNibydz(8);   ...             
       dNibydx(1)   dNibydy(1)           0             dNibydx(2)   dNibydy(2)          0               dNibydx(3)    dNibydy(3)            0             dNibydx(4)    dNibydy(4)        0            ...
       dNibydx(5)   dNibydy(5)           0             dNibydx(6)   dNibydy(6)          0               dNibydx(7)    dNibydy(7)            0             dNibydx(8)    dNibydy(8)        0];
  
BTDBJ1=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{ -1/sqrt(3), -1/sqrt(3), -1/sqrt(3)});
BTDBJ2=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{  1/sqrt(3), -1/sqrt(3), -1/sqrt(3)});
BTDBJ3=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{  1/sqrt(3),  1/sqrt(3), -1/sqrt(3)});
BTDBJ4=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{ -1/sqrt(3),  1/sqrt(3), -1/sqrt(3)});

BTDBJ5=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{ -1/sqrt(3), -1/sqrt(3), 1/sqrt(3)});
BTDBJ6=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{  1/sqrt(3), -1/sqrt(3), 1/sqrt(3)});
BTDBJ7=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{  1/sqrt(3),  1/sqrt(3), 1/sqrt(3)});
BTDBJ8=subs(B'*D*B*det(J),{Xi, Eta, Zeta},{ -1/sqrt(3),  1/sqrt(3), 1/sqrt(3)});



k=BTDBJ1+BTDBJ2+BTDBJ3+BTDBJ4+BTDBJ5+BTDBJ6+BTDBJ7+BTDBJ8;





   