matlab四节点矩形刚度矩阵代码

%该程序计算平面四节点四边形等参元的刚度矩阵 %输入弹性模量E、泊松比NU和厚度h %输入四个节点1、2、3、4的坐标x1,x2,x3,x3,y1,y2,y3,y4 %输入平面问题性质参数ID（1为平面应力，2为平面应变） %输出单元刚度矩阵

%%%%%%%%%%%%%%%%%%%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