matlab有限元需要的几个基本函数

五个函数如下:

function [L,theta] = PTE_length(x1,y1,x2,y2)
%根据给出的坐标（x1.y1）,(x2,y2)计算返回单元的长度与角度
L=sqrt((x2-x1)^2+(y2-y1)^2);
theta=atan((y2-y1)/(x2-x1));
end

function [k] = PTE_stiffness(E,A,L,theta)
%根据给出的单元弹性模量E，横街面积A， 长度L和角度theta计算返回4*4单元刚度矩阵k
T=[cos(theta) sin(theta) 0 0;0 0 cos(theta) sin(theta)];
k=[1 -1;-1 1];
k=k*E*A/L;
k=T'*k*T;
end

function [stress] = PTE_stress(E,L,theta,q)
%根据给出的单元弹性模量E，横街面积A，长度L和角度theta以及整体坐标系下单元节点位移向量u（4*1）计算单元应力；
T=[cos(theta) sin(theta) 0 0;0 0 cos(theta) sin(theta)];
stress=E/L*[-1 1]*T*q;
end

function [P] = PTE_force(E,A,L,theta,q)
%根据给出的单元弹性模量E，横街面积A，长度L，角度theta以及整体坐标系下单元节点位移向量u（4*1）,计算返回整体坐标系下的单元节点力P
P=PTE_stiffness(E,A,L,theta)*q;
End

function [K] = PT_assemble(K,k,i,j)
%将连接节点i，j的单元门刚度矩阵k集成到整体刚度矩阵K。计算返回2n*2n整体刚度矩阵K
K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);
K(2*i-1,2*i)  =K(2*i-1,2*i)  +k(1,2);
K(2*i,2*i-1)  =K(2*i,2*i-1)  +k(2,1);
K(2*i,2*i)    =K(2*i,2*i)    +k(2,2);

K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);
K(2*i-1,2*j)  =K(2*i-1,2*j)  +k(1,4);
K(2*i,2*j-1)  =K(2*i,2*j-1)  +k(2,3);
K(2*i,2*j)    =K(2*i,2*j)    +k(2,4);

K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);
K(2*j-1,2*i)  =K(2*j-1,2*i)  +k(3,2);
K(2*j,2*i-1)  =K(2*j,2*i-1)  +k(4,1);
K(2*j,2*i)    =K(2*j,2*i)    +k(4,2);

K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);
K(2*j-1,2*j)  =K(2*j-1,2*j)  +k(3,4);
K(2*j,2*j-1)  =K(2*j,2*j-1)  +k(4,3);
K(2*j,2*j)    =K(2*j,2*j)    +k(4,4);
End