clc,clear,close all
tic
nodes=zeros(3,244); %创建节点
ele0=zeros(2,5);ele1=zeros(24,5);ele2=zeros(60,5);ele3=zeros(20,5);
ele4=zeros(64,5);ele5=zeros(168,5);ele6=zeros(28,5);ele7=zeros(144,5);
ele8=zeros(396,5);ele9=zeros(36,5);
d1=4270/2;d2=5335;r2=4265;d3=7470;r3=6400;d4=9605;r4=8535;
l3=43890;l2=l3+29260;l1=l2+21950;
dl3=43890/12;dl2=29260/8;dl1=21950/6;
coordinate2=[d2,r2,0,-r2,-d2,-r2,0,r2;0,r2,d2,r2,0,-r2,-d2,-r2];%定义坐标
coordinate3=[d3,r3,d1,-d1,-r3,-d3,-d3,-r3,-d1,d1,r3,d3;...
d1,r3,d3,d3,r3,d1,-d1,-r3,-d3,-d3,-r3,-d1];
coordinate4=[d4,r4,d1,-d1,-r4,-d4,-d4,-r4,-d1,d1,r4,d4;...
d1,r4,d4,d4,r4,d1,-d1,-r4,-d4,-d4,-r4,-d1];
%节点坐标
for i=1:6
    nodes(:,4*i-3:4*i)=[d1,-d1,-d1,d1;d1,d1,-d1,-d1;l1-(i-1)*dl1*ones(1,4)];
end
for i=1:8
     nodes(:,24+8*(i-1)+1:24+8*i)=[coordinate2;l2-(i-1)*dl2*ones(1,8)];
end
for i=1:12
    nodes(:,88+12*(i-1)+1:88+12*i)=[coordinate3;l3-(i-1)*dl3*ones(1,12)];
end
nodes(:,233:244)=[coordinate4;zeros(1,12)];
%单元信息：单元编号，单元节点，单元节点
ele0(:,2:3)=[1,3;2,4];
for i=1:6
   ele1(4*i-3:4*i,2:3)=[4*i-3,4*i-2;4*i-2,4*i-1;4*i-1,4*i;4*i,4*i-3];%1-24
   if i==6
      ele3(1:5,2:3)=[21,32;21,25;21,26;21,27;21,28];
      ele3(6:10,2:3)=[22,26;22,27;22,28;22,29;22,30];
      ele3(11:15,2:3)=[23,28;23,29;23,30;23,31;23,32];
      ele3(16:20,2:3)=[24,30;24,31;24,32;24,25;24,26];
   else
      ele2(12*(i-1)+3-2:12*(i-1)+3,2:3)=...
[4*(i-1)+1,4*i+4; 4*(i-1)+1,4*i+1; 4*(i-1)+1,4*i+2];
      for j=2:3
         ele2(12*(i-1)+3*j-2:12*(i-1)+3*j,2:3)=...
[4*(i-1)+j,4*i+j-1; 4*(i-1)+j,4*i+j; 4*(i-1)+j,4*i+j+1];
      end
      ele2(12*(i-1)+12-2:12*(i-1)+12,2:3)=[4*(i-1)+4,4*i+4-1;...
4*(i-1)+4,4*i+4; 4*(i-1)+4,4*i+1];
   end
end
for i=1:8
   ele4(8*i-7:8*i,2:3)=[8*i-7,8*i-6;8*i-6,8*i-5;8*i-5,8*i-4;...
8*i-4,8*i-3;8*i-3,8*i-2;8*i-2,8*i-1;8*i-1,8*i;8*i,8*i-7;];
   if i==8
      ele6(1:4,2:3)=[81,99;81,100;81,89;81,90]; 
      for j=1:3
         ele6(4*(j-1)+5:4*(j-1)+8,2:3)=[81+2*j,88+3*j-1;...
81+2*j,88+3*j;81+2*j,88+3*j+1;81+2*j,88+3*j+2]; 
      end  
      for j=1:4
         ele6(3*(j-1)+17:3*(j-1)+19,2:3)=[80+2*j,89+3*(j-1);...
80+2*j,89+3*(j-1)+1;80+2*j,89+3*(j-1)+2]; 
      end  
   else
      ele5(24*(i-1)+1:24*(i-1)+3,2:3)=...
[8*(i-1)+1,8*i+8;8*(i-1)+1,8*i+1; 8*(i-1)+1,8*i+2];
      for j=2:7
          ele5(24*(i-1)+3*j-2:24*(i-1)+3*j,2:3)=...
[8*(i-1)+j,8*i+j-1; 8*(i-1)+j,8*i+j; 8*(i-1)+j,8*i+j+1];
      end
      ele5(24*(i-1)+24-2:24*(i-1)+24,2:3)=...
[8*(i-1)+8,8*i+8-1; 8*(i-1)+8,8*i+8; 8*(i-1)+8,8*i+1];
   end
end
ele4=ele4+24;
ele5=ele5+24;
for i=1:12
   ele7(12*i-11:12*i,2:3)=[12*i-11,12*i-10;12*i-10,12*i-9;12*i-9,...
12*i-8;12*i-8,12*i-7;12*i-7,12*i-6;12*i-6,12*i-5;12*i-5,12*i-4;...
12*i-4,12*i-3;12*i-3,12*i-2;12*i-2,12*i-1;12*i-1,12*i;12*i,12*i-11;];%1-24
   if i==12
      ele9(1:3,2:3)=[233,232;233,221;233,222]; 
      for j=2:11
         ele9(3*(j-1)+1:3*(j-1)+3,2:3)=[232+j,220+j-1;232+j,220+j;232+j,220+j+1]; 
      end  
      ele9(34:36,2:3)=[244,231;244,232;244,221]; 
   else
      ele8(36*(i-1)+1:36*(i-1)+3,2:3)=[12*(i-1)+1,...
12*i+12;12*(i-1)+1,12*i+1; 12*(i-1)+1,12*i+2];
      for j=2:11
          ele8(36*(i-1)+3*(j-1)+1:36*(i-1)+3*j,2:3)=[12*(i-1)+j,...
12*i+j-1; 12*(i-1)+j,12*i+j; 12*(i-1)+j,12*i+j+1];
      end
      ele8(36*(i-1)+36-2:36*(i-1)+36,2:3)=[12*(i-1)+12,12*i+12-1; ...
12*(i-1)+12,12*i+12; 12*(i-1)+12,12*i+1];
   end
end
ele7=ele7+88;
ele8=ele8+88;
ele=[ele0;ele1;ele2;ele3;ele4;ele5;ele6;ele7;ele8;ele9];
for i=1:length(ele(:,1))
   ele(i,1)=i;
end
A=4;E=2.1E005; % 定义各单元截面面积和材料弹性模量
ele(:,4:5)=ones(length(ele(:,1)),2)*[A 0;0 E]; %定义单元信息: 截面积，弹性模量
Load=[1 0 400 -100;2 0 400 -100]; %载荷信息：单元编号，x方向力，y方向力，z方向力
x=nodes(1,:);y=nodes(2,:);z=nodes(3,:);
Constr=zeros(12,4);%定义约束
for i=1:12
   Constr(i,1)=i+232;
end
[U, Strain, Stress, AxialForce]=Truss3DStaticsFEA(x,y,z,ele,Load,Constr,10) 
toc