main.m

基于力密度法的索杆张力结构找形程序

function main
%以坐标的增量，有限元程序


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  PART  1    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%开始，读入数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nodenum,nodefree,nodefix,cab_num,node_all,nfix_n,cab_all,node_no,node_x,node_y,node_z,px,py,pz,cab_n,cab_i,cab_j,E,Real,Area,Etype,barinfor,filename1]=fiinput;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%读入杆件初始长度s0%%%%%%%%%%%%%%%%%%%%%%%%5%%%
[s0,filename2]=intiallength(cab_num);

% nnn=24;         %等分数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  PART  2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%拓扑矩阵分离自由拓扑C 固定节点的拓扑 Cf 以及自由节点的节点力%%%%%%%%%%%%%%%%%
[C1,C,Cf,pxl,pyl,pzl,bianhao]=topu(cab_num,cab_n,cab_i,cab_j,nfix_n,nodefix,nodenum,node_no,px,py,pz);
%%%%%%%%%分离自由节点 l 和约束节点 fl 以及自由和约束节点的节点号%%%%%%%%%%%%%%%%%%%%%%%
[Xf,X,Yf,Y,Zf,Z,node_fixedno,node_freeno]=nodefenli(nodefix,nodenum,node_no,nfix_n,node_x,node_y,node_z);

%%%%%%%%%%%%%%%%%%%%%%%%确定初始迭代的坐标x0,y0,z0只需要给定自由节点的坐标即可%%%%%%%%%%%%%%%%%%%%%%%%%%%%
weight=7850*9.8*Area.*s0;

% % 单索的形状方程
% syms f l x c t E A r w s0 ux uz L
% z=4*f*x*(l-x)/l^2+c*x/l;
% z=diff(z,x);
% SS=int(sqrt(z^2+1),x,0,l);
% derts=t*l/(r*E*A)*int(1+z^2,x,0,l)
% derts=subs(derts,f,w*r/(8*t))
% simple(derts)
% SS=subs(SS,f,w*r/(8*t));
% SS=subs(SS,l,L);
% SS=simple(SS);
% SS=simple(SS);
% SS=simple(SS)
% [y,sigama]=subexpr(SS)
% FS=SS-s0-t/(E*A)*(l^2/r+c^2/r+w^2*r/12/t^2)
% 
% 
% 
% 
% FSpianux=diff(FS,r)*ux/r+diff(FS,l)*ux/l
% 
% 
% 
% FSpianuz=diff(FS,r)*uz/r+diff(FS,l)*uz/l+diff(FS,c)*sign(uz)
% 
% 
% 
% FSpiant=diff(FS,t)

% weight=0.0001*s0
%迭代步长
stepp=0.5;
time=0;
% xx0=X;yy0=Y;zz0=Z;xx1=X+10;yy1=Y+10;zz1=Z+10;xx2=X+100;yy2=Y+100;zz2=Z+100;xx3=X+50;yy3=Y+50;zz3=Z+50;
%%%%%%%%%%%%%%%%%%%%%%%%%%计算迭带的参数fx(x,y,z,h),以及一阶增量%%%%%%%%%%%
while 1>0
    ux=C*X+Cf*Xf;
    uy=C*Y+Cf*Yf;
    uz=C*Z+Cf*Zf;
    l=(ux.^2+uy.^2).^0.5;
    r=(ux.^2+uy.^2+uz.^2).^0.5;
    for i=1:1:cab_num
        L(i,i)=l(i);
        Ux(i,i)=ux(i);
        Uy(i,i)=uy(i);
        Uz(i,i)=uz(i);
        R(i,i)=r(i);
    end
%   [cab_all(:,1),l,r,uz,s0]
    %%%%%%%%%%%%%%%%%%%%%%%确定在初始态的索的水平分力的H的大小%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%区别压杆和拉索的t力不同%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for i=1:1:cab_num
        if Etype(i)==0
            tforce(i,1)=Tbar(r(i),s0(i),E(i),Area(i));
        elseif Etype(i)==1
            tforce(i,1)=Tcable(ux(i),uy(i),uz(i),s0(i),E(i),Area(i),weight(i));
        else
            errordlg({'单元类型有错误，需检查'});
        end
    end

%
%   [cab_all(:,1),l,r,uz,s0,tforce/1000000,E/1E11,Area];
%
    for i=1:1:cab_num
        Tforce(i,i)=tforce(i,1);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%计算变量fx0,fy0,fz0和各个gama系数,继而算出坐标增量%%%%%%%%%
    [fx0,fy0,fz0]=ffct(C,Ux,Uy,Uz,R,tforce);
    [gamarx,gamary,gamarz]=gamar_fct(R,L,Ux,Uz,Uy,Tforce,E,Area,weight,s0,cab_num,Etype);
    if rank(C'*gamarx*C)~=nodefree|rank(C'*gamary*C)~=nodefree|rank(C'*gamarz*C)~=nodefree
        errordlg({'gama的为非满秩矩阵，需要进行检验'});
        break;
    end
    dertX=inv(C'*gamarx*C)*(pxl-fx0)
    dertY=inv(C'*gamary*C)*(pyl-fy0)
    dertZ=inv(C'*gamarz*C)*(pzl-0.5*abs(C')*weight-fz0)
    if max(abs(dertX))<0.000005&max(abs(dertY))<0.000005&max(abs(dertZ))<0.000005
        X=X+dertX;
        Y=Y+dertY;
        Z=Z+dertZ;
        break
    end
    steppx=stepp/1;steppy=stepp/1;
    for i=1:1:nodefree
        if dertX(i)>steppx
            dertX(i)=steppx;
        elseif dertX(i)<-steppx
            dertX(i)=-steppx;
        end
        if dertY(i)>steppy
            dertY(i)=steppy;
        elseif dertY(i)<-steppy
            dertY(i)=-steppy;
        end
        if dertZ(i)>stepp
            dertZ(i)=stepp;
        elseif dertZ(i)<-stepp
            dertZ(i)=-stepp;
        end
    end
    X=X+dertX;
    Y=Y+dertY;
    Z=Z+dertZ
    if mod(time,4)==0
        xx0=X;
        yy0=Y;
        zz0=Z;
    elseif mod(time,4)==1
        xx1=X;
        yy1=Y;
        zz1=Z;
    elseif mod(time,4)==2
        xx2=X;
        yy2=Y;
        zz2=Z;
    elseif mod(time,4)==3
        xx3=X;
        yy3=Y;
        zz3=Z;
    end
    if mod(time,4)==3
%         fprintf('_____________________________')
        stepp=Step(xx0,yy0,zz0,xx1,yy1,zz1,xx2,yy2,zz2,xx3,yy3,zz3,stepp);
    end
    time=time+1
    stepp
%     if mod(time,40)==0
%         ux=C*X+Cf*Xf;
%         uy=C*Y+Cf*Yf;
%         uz=C*Z+Cf*Zf;
%         l=(ux.^2+uy.^2).^0.5;
%         r=(ux.^2+uy.^2+uz.^2).^0.5;
%         for i=1:1:cab_num
%             if Etype(i)==0
%                 tforce(i,1)=Tbar(r(i),s0(i),E(i),Area(i));
%             elseif Etype(i)==1
%                 tforce(i,1)=Tcable(ux(i),uy(i),uz(i),s0(i),E(i),Area(i),weight(i));
%             else
%                 errordlg({'单元类型有错误，需检查'});
%             end
%         end
%         Xall=[X',Xf']';
%         Yall=[Y',Yf']';
%         Zall=[Z',Zf']';
%         outputs(tforce,E,Area,l,uz,weight,cab_num,Xall,Yall,Zall,r,s0,Etype,bianhao,cab_n,filename1)
%         pictoautocad(Xall,Yall,Zall,C1,weight,l,tforce,Etype,r,filename2)
%         todate(node_all,nodenum,nodefree,nodefix,cab_num,Xall,Yall,Zall,bianhao,nfix_n,cab_all,filename1)
%     end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%输出长度s和变换长度s0%%%%%%%%%%%%%%%%%%%%%
ux=C*X+Cf*Xf;
uy=C*Y+Cf*Yf;
uz=C*Z+Cf*Zf;
l=(ux.^2+uy.^2).^0.5;
r=(ux.^2+uy.^2+uz.^2).^0.5;
for i=1:1:cab_num
    if Etype(i)==0
        tforce(i,1)=Tbar(r(i),s0(i),E(i),Area(i));
    elseif Etype(i)==1
        tforce(i,1)=Tcable(ux(i),uy(i),uz(i),s0(i),E(i),Area(i),weight(i));
    else
        errordlg({'单元类型有错误，需检查'});
    end
end
Xall=[X',Xf']';
Yall=[Y',Yf']';
Zall=[Z',Zf']'
outputs(tforce,E,Area,l,uz,weight,cab_num,Xall,Yall,Zall,r,s0,Etype,bianhao,cab_n,filename1)
pictoautocad(Xall,Yall,Zall,C1,weight,l,tforce,Etype,r,filename2)
todate(node_all,nodenum,nodefree,nodefix,cab_num,Xall,Yall,Zall,bianhao,nfix_n,cab_all,filename1)
% X
% Y
% Z
% Hforce
% weight