function w=brickelementstiffness(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;
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;
J31=diff(X,r);J32=diff(Y,r);J33=diff(Z,r);
J11=diff(X,s);J12=diff(Y,s);J13=diff(Z,s);
J21=diff(X,t);J22=diff(Y,t);J23=diff(Z,t);
J=det([J11 J12 J13;J21 J22 J23;J31 J32 J33]);

difX=det([diff(N1,s) J12 J13;diff(N1,t) J22 J23;diff(N1,r) J32 J33])/J;
difY=det([J11 diff(N1,s) J13;J21 diff(N1,r) J23;J31 diff(N1,r) J33])/J;
difZ=det([J11 J12 diff(N1,s);J21 J22 diff(N1,r);J31 J32 diff(N1,r)])/J;
B1=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N2,s) J12 J13;diff(N2,t) J22 J23;diff(N2,r) J32 J33])/J;
difY=det([J11 diff(N2,s) J13;J21 diff(N2,r) J23;J31 diff(N2,r) J33])/J;
difZ=det([J11 J12 diff(N2,s);J21 J22 diff(N2,r);J31 J32 diff(N2,r)])/J;
B2=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N3,s) J12 J13;diff(N3,t) J22 J23;diff(N3,r) J32 J33])/J;
difY=det([J11 diff(N3,s) J13;J21 diff(N3,r) J23;J31 diff(N3,r) J33])/J;
difZ=det([J11 J12 diff(N3,s);J21 J22 diff(N3,r);J31 J32 diff(N3,r)])/J;
B3=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N4,s) J12 J13;diff(N4,t) J22 J23;diff(N4,r) J32 J33])/J;
difY=det([J11 diff(N4,s) J13;J21 diff(N4,r) J23;J31 diff(N4,r) J33])/J;
difZ=det([J11 J12 diff(N4,s);J21 J22 diff(N4,r);J31 J32 diff(N4,r)])/J;
B4=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N5,s) J12 J13;diff(N5,t) J22 J23;diff(N5,r) J32 J33])/J;
difY=det([J11 diff(N5,s) J13;J21 diff(N5,r) J23;J31 diff(N5,r) J33])/J;
difZ=det([J11 J12 diff(N5,s);J21 J22 diff(N5,r);J31 J32 diff(N5,r)])/J;
B5=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N6,s) J12 J13;diff(N6,t) J22 J23;diff(N6,r) J32 J33])/J;
difY=det([J11 diff(N6,s) J13;J21 diff(N6,r) J23;J31 diff(N6,r) J33])/J;
difZ=det([J11 J12 diff(N6,s);J21 J22 diff(N6,r);J31 J32 diff(N6,r)])/J;
B6=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N7,s) J12 J13;diff(N7,t) J22 J23;diff(N7,r) J32 J33])/J;
difY=det([J11 diff(N7,s) J13;J21 diff(N7,r) J23;J31 diff(N7,r) J33])/J;
difZ=det([J11 J12 diff(N7,s);J21 J22 diff(N7,r);J31 J32 diff(N7,r)])/J;
B7=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];

difX=det([diff(N8,s) J12 J13;diff(N8,t) J22 J23;diff(N8,r) J32 J33])/J;
difY=det([J11 diff(N8,s) J13;J21 diff(N8,r) J23;J31 diff(N8,r) J33])/J;
difZ=det([J11 J12 diff(N8,s);J21 J22 diff(N8,r);J31 J32 diff(N8,r)])/J;
B8=[difX 0 0; 0 difY 0; 0 0 difZ;difX difY 0;difX 0 difZ;0 difY difZ];
B=[B1 B2 B3 B4 B5 B6 B7 B8];
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];
BD=transpose(B)*D*B*J;
K=int(int(int(BD,t,-1,1),r,-1,1),s,-1,1);
w=double(K);