clc,clear
tic
d1=37.5;d2=37.5;d3=100;h1=200;h2=100;
x=[-d1 d1 -d2 d2 d2 -d2 -d3 d3  d3 -d3]; % 节点x轴方向坐标
y=[0 0 d2 d2 -d2 -d2 h2 d3 -d3 -d3]; % 节点y轴方向坐标
z=[h1 h1 h2 h2 h2 h2 0 0 0 0]; % 节点z轴方向坐标
A=3;E=2.1E005; % 定义截面面积和弹性模量
%单元信息：编号，节点1编号，节点2编号，截面面积，弹性模量
ele=[1 1 2 A E;2, 4, 1 A E;3 2 3 A E;4, 2, 6 A E;5, 1, 5 A E;6, 2, 4 A E;7, 2, 5 A E;...
   8, 1, 3 A E;9, 1, 6 A E;10, 3, 6 A E;11, 5, 4 A E;12, 3, 4 A E;13, 5, 6 A E;...
   14, 3, 10 A E;15, 7, 6 A E;16, 9, 4 A E;17, 8, 5 A E;18, 7, 4 A E;19, 3, 8 A E;...
   20, 10, 5 A E;21, 9, 6 A E;22, 10, 6 A E;23, 3, 7 A E;24, 8, 4 A E;25, 9, 5 A E];
%载荷信息：节点编号，x方向力，y方向力，z方向力
Load=[1 0 200 0;2 0 200 0];
%约束：节点编号，x方向约束，y方向约束，z方向约束（未约束填NaN）
Constr=[7 0 0 0;8 0 0 0;9 0 0 0;10 0 0 0];
[U, Strain, Stress, AxialForce]=Truss3DStaticsFEA(x,y,z,ele,Load,Constr,40) 
toc