%此程序计算杆的总刚度矩阵及节点位移，分五步完成：输入各单元数据、计算单元刚度矩阵

%组集总刚度矩阵、计算输出总刚度矩阵 、计算输出节点位移

%2011.4

%输入个单元数据

%输入单元节点编号矩阵，每一行代表该单元的节点编号，即用分号将各单元分开

%用逗号将每个单元内的节点分开，注意输入节点顺序要相同，即各单元要么同为顺时针

%，要么同为逆时针，按单元编号顺序排列

%如[1,2,3;4,3,1]表示两个三角形中的节点编号分别为(1,2,3)、(4,3,1)

cod=input('please input the node of each element in order:');

%计算单元个数，nm为单元个数

[nm,nmn]=size(cod);

%输入节点坐标，每一行代表该节点的坐标，按节点编号顺序排列，即用分号将节点分开

%用逗号将每个节点的坐标分开，按单元编号顺序排列

%如[1,2;2,3;1,3]表示三个节点的横纵坐标分别为(1,2)、(2,3)、(1,3)

con=input('please input the coordinates (m) of each node in order:');

%计算结点个数，nn为结点个数

[nn,nnn]=size(con);

%输入单元的弹性模量

E=1e9*input('please input E (GPa) of each element in order:');

%输入泊松比miu

miu=input('please input the poison ratio miu of the material:');

%输入厚度t

t=input('please input thickness t (m) of the material:');

%输入节点载荷及相应的节点坐标

P=input('please input the force (N) on the node(if the force is unknown please input nan):');

%输入节点位移

Displacement=input('please input the displacement (m) of the node(if the displacement is unknown please input nan):');

%输入解决问题的性质，平面应力问题性质输入1，平面应变问题输入2

disp('please input the character of the problem:')

ID=input('For Plane stress problem input 1,for Plane strain problem input 2:');

%计算单元刚度矩阵

kele=zeros(6,6,nm);

for i=1:nm

x1=con(cod(i,1),1);

y1=con(cod(i,1),2);

x2=con(cod(i,2),1);

y2=con(cod(i,2),2);

x3=con(cod(i,3),1);

y3=con(cod(i,3),2);

kele(:,:,i)=Triangle_Stiffness(E,miu,t,x1,y1,x2,y2,x3,y3,ID);%ID输入解决问题的性质

end

%组集总刚度矩阵

%定义空的总刚矩阵

kg=sparse(2*nn,2*nn);%定定义稀疏矩阵，只储存非零元素

%组集各单元矩阵

for n=1:nm

num1=cod(n,1);%提取节点坐标

num2=cod(n,2);

num3=cod(n,3);

kg=Triangle_Assembly(kg,kele(:,:,n),num1,num2,num3);%此函数中实现变半带宽储存

end

%输出总刚度矩阵

disp('kg=')

disp(kg)

%计算节点位移

kg=kg+tril(kg,-1)';%把下三角矩阵变为整体矩阵

disp('kg_full=')

disp(full(kg))

%置一赋零法引入边界条