truss.m

6根杆四个节点

         er=1.;
         ea=6e-4;
         ee=210e9; 
         ne=2;
         nfe=4;
         nf=2;
         p=[1.e6;1.e6];
         xy=[0,0,0.707,0.707;0,1.414,0.707,0.707];
         ndg=[3,3,1,2;3,3,1,2];
         sk=zeros(nf,nf);
         sm=zeros(nf,nf);
         m=zeros(4);
         for kk=1:ne,
         ek=zeros(4,4);
         em=zeros(4,4);
         xkk1=xy(kk,1);
         ykk1=xy(kk,2);
         xkk2=xy(kk,3);
         ykk2=xy(kk,4);
         lkk=(xkk1-xkk2)^2+(ykk1-ykk2)^2;
         lkk=sqrt(lkk);
         ca=(xkk2-xkk1)/lkk;
         sa=(ykk2-ykk1)/lkk;
         ca2=ca*ca;
         sa2=sa*sa;
         csa=ca*sa;
         ck=ee*ea/lkk;
         cm=er*ea*lkk/2.;
         ek=[ca2,csa,-ca2,-csa;csa,sa2,-csa,-sa2;-ca2,-csa,ca2,csa;-csa,-sa2,csa,sa2];
         ek=ck*ek;
%         em=[ca2,csa,0,0;csa,sa2,0,0;0,0,ca2,csa;0,0,csa,sa2];
         em=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];
         em=cm*em;
%        p=[ca,-sa,0,0;sa,ca,0,0;0,0,ca,-sa;0,0,sa,ca];
%        ek=p*ek*p';
%        em=p*em*p';
            for i=1:nfe,
                 m(i)=ndg(kk,i);
             end
             for i=1:nfe,
                 mi=m(i);
                 if(mi<=nf),
                 for j=1:nfe,
                     mj=m(j);
                      if(mj<=nf),
                        sk(mi,mj)=sk(mi,mj)+ek(i,j);
                        sm(mi,mj)=sm(mi,mj)+em(i,j);
                      end
                   end
                  end
             end
         end
uv=inv(sk)*p
%         [v,d]=eig(sk,sm);
    %     ev=diag(d);
       %  ev=sqrt(ev)
       %  [evs,ser]=sort(ev)
       %  v
        % for j=1:nf,
        %     for i=1:nf,
       %      vec(i,j)=v(i,ser(j));
         %    end
        % end
        % vec