%-----------------------------------------------------------
function?k=Hexahedral3D8Node_Stiffness（ENUx1y1z1x2y2z2x3y3z3x4y4z4x5y5z5x6y6z6x7y7z7x8y8z8）
%該函數計算單元的剛度矩陣
%輸入彈性模量E，泊松比NU
%輸入8個節點的坐標x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4、
%輸入x5、y5、z5、x6、y6、z6、x7、y7、z7、x8、y8、z8
%輸出單元剛度矩陣k（24X24）?
%---------------------------------------------------------------
syms?s?t?n;??????????????????%?定義局部坐標系?
%???定義形狀函數?N
N1=（1+s）*（1-t）*（1-n）/8;
N2=（1+s）*（1+t）*（1-n）/8;
N3=（1-s）*（1+t）*（1-n）/8;
N4=（1-s）*（1-t）*（1-n）/8;
N5=（1+s）*（1-t）*（1+n）/8;
N6=（1+s）*（1+t）*（1+n）/8;
N7=（1-s）*（1+t）*（1+n）/8;
N8=（1-s）*（1-t）*（1+n）/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;
%???定義雅可比矩陣?
J=[diff（xs）diff（ys）diff（zs）;diff（xt）diff（yt）diff（zt）;diff（xn）diff（yn）diff（zn）];
Jdet=det（J）;
J;
Jdet;
%???定義?B?矩陣的系數
a=diff（yt）*diff（zn）-diff（zt）*diff（yn）;
b=diff（ys）*diff（zn）-diff（zs）*diff（yn）;
c=diff（ys）*diff（zt）-diff（zs）*diff（yt）;
d=diff（xt）*diff（zn）-diff（zt）*diff（xn）;
e=diff（xs）*diff（zn）-diff（zs）*diff（xn