%%%%%%%%%%%%%%%%%%%Quad2D4NodeGuass_Stiffness%%%%begin%%%%%%%%%%
function?k=Quad2D4NodeGuass_Stiffness（ENUhx1y1x2y2x3y3x4y4ID）
%該程序計算平面四節點四邊形等參元的剛度矩陣
%輸入彈性模量E、泊松比NU和厚度h
%輸入四個節點1、2、3、4的坐標x1x2x3x3y1y2y3y4
%輸入平面問題性質參數ID（1為平面應力，2為平面應變）
%輸出單元剛度矩陣
%-------------------------------------------------------------------
Ai=1
Aj=1;

syms?s?t;
a=[-（1-t）*x1+（1-t）*x2+（1+t）*x3-（1+t）*x4]/4;
b=[-（1-s）*x1-（1+s）*x2+（1+s）*x3+（1-s）*x4]/4;
c=[-（1-t）*y1+（1-t）*y2+（1+t）*y3-（1+t）*y4]/4;
d=[-（1-s）*y1-（1+s）*y2+（1+s）*y3+（1-s）*y4]/4;
J=[a?c;b?d];
L=det（J）;
H=[d