function y=Brick8ElementStiffness1(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8)
syms r s t;% r-x s-y t-z
N1=(1+r)*(1-s)*(1-t)/8;
N2=(1+r)*(1+s)*(1-t)/8;
N3=(1-r)*(1+s)*(1-t)/8;
N4=(1-r)*(1-s)*(1-t)/8;
N5=(1+r)*(1-s)*(1+t)/8;
N6=(1+r)*(1+s)*(1+t)/8;
N7=(1-r)*(1+s)*(1+t)/8;
N8=(1-r)*(1-s)*(1+t)/8;
X=N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8;
Y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8;
Z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8;
J11=diff(X,r);J12=diff(Y,r);J13=diff(Z,r);
J21=diff(X,s);J22=diff(Y,s);J23=diff(Z,s);
J31=diff(X,t);J32=diff(Y,t);J33=diff(Z,t);
J=[J11 J12 J13;J21 J22 J23;J31 J32 J33];
N=[N1 N2 N3 N4 N5 N6 N7 N8];B=zeros(6,24);B=sym(B);
for i=1:8
difN=inv(J)*[diff(N(i),r);diff(N(i),s);diff(N(i),t)];
B(1:6,(3*i-2):3*i)=[difN(1) 0 0; 0 difN(2) 0; 0 0 difN(3);difN(2) difN(1) 0;...
    0 difN(3) difN(2);difN(3) 0 difN(1)];
end
D=E/((1-2*NU)*(1+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];
J=det(J);
BD=transpose(B)*D*B*J;
K=int(int(int(BD,t,-1,1),r,-1,1),s,-1,1);
y=double(K);