k_assembly.m

FEM tools for caculation of nonlinear problems

        if in_data.mater.E~=0    E1=in_data.mater.E;      else E1=in_data.EL(i,5);    end;
        if in_data.mater.G~=0    G1=in_data.mater.G;      else G1 = in_data.EL(i,6);  end;
        if in_data.mater.b~=0    b1=in_data.mater.b;      else b1 = in_data.EL(i,7);  end;
        if in_data.mater.ho~=0   ho1=in_data.mater.ho;    else ho1 = in_data.EL(i,8); end;
        if in_data.mater.rho~=0  rho_=in_data.mater.rho;  else rho_= in_data.EL(i,9); end;

        [M_loc,Klc] = D3_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node1,4), ...
              in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node2,4),E1,G1,b1,ho1,rho_);

        t1 = node1*dofPnd; t2 = node2*dofPnd;

        bg = [(t1-(dofPnd-1)):t1]; en = [(t2-(dofPnd-1)):t2];
        
        Kgl([bg en],[bg en]) = Kgl([bg en],[bg en]) + Klc;
        M([bg en],[bg en])   = M([bg en],[bg en])   + M_loc;
    end;
   
     if in_data.EL(i,2)==6
         
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
        node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
        node5 = find(in_data.ND(:,1)==in_data.EL(i,7));
        node6 = find(in_data.ND(:,1)==in_data.EL(i,8));
        node7 = find(in_data.ND(:,1)==in_data.EL(i,9));
        node8 = find(in_data.ND(:,1)==in_data.EL(i,10));
        
        x1 = in_data.ND(node1,2); y1 = in_data.ND(node1,3); z1 = in_data.ND(node1,4);
        x2 = in_data.ND(node2,2); y2 = in_data.ND(node2,3); z2 = in_data.ND(node2,4);
        x3 = in_data.ND(node3,2); y3 = in_data.ND(node3,3); z3 = in_data.ND(node3,4);
        x4 = in_data.ND(node4,2); y4 = in_data.ND(node4,3); z4 = in_data.ND(node4,4);
        x5 = in_data.ND(node5,2); y5 = in_data.ND(node5,3); z5 = in_data.ND(node5,4);
        x6 = in_data.ND(node6,2); y6 = in_data.ND(node6,3); z6 = in_data.ND(node6,4);
        x7 = in_data.ND(node7,2); y7 = in_data.ND(node7,3); z7 = in_data.ND(node7,4);
        x8 = in_data.ND(node8,2); y8 = in_data.ND(node8,3); z8 = in_data.ND(node8,4);

        if in_data.mater.E~=0    Em=in_data.mater.E;      else  Em=in_data.EL(i,11);       end;
        if in_data.mater.miu~=0  miu_1=in_data.mater.miu; else  miu_1 = in_data.EL(i,12);  end;
        if in_data.mater.rho~=0  rho_=in_data.mater.rho;  else  rho_=in_data.EL(i,13);     end;
        
        [Klc,Bsys,Esys,V] = D3_LB (x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,...
            x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,Em,miu_1);

        sf = ones(24,1);
        fM = V*rho_/8.*sf';

        t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd; t4 = node4*dofPnd; 
        t5 = node5*dofPnd; t6 = node6*dofPnd; t7 = node7*dofPnd; t8 = node8*dofPnd;

        bg = [(t1-(dofPnd-1)):t1];
        md1= [(t2-(dofPnd-1)):t2];
        md2= [(t3-(dofPnd-1)):t3];
        md3= [(t4-(dofPnd-1)):t4];
        md4= [(t5-(dofPnd-1)):t5];
        md5= [(t6-(dofPnd-1)):t6];
        md6= [(t7-(dofPnd-1)):t7];
        en = [(t8-(dofPnd-1)):t8];
        
        Kgl([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) = ...
            Kgl([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) + Klc;
        M([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) = ...
            M([bg md1 md2 md3 md4 md5 md6 en],[bg md1 md2 md3 md4 md5 md6 en]) + diag(fM);
    end;

     if in_data.EL(i,2)==10
        
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
        node4 = find(in_data.ND(:,1)==in_data.EL(i,6));

        if in_data.mater.E~=0    Em=in_data.mater.E;      else  Em=in_data.EL(i,7);       end;
        if in_data.mater.miu~=0  miu_1=in_data.mater.miu; else  miu_1 = in_data.EL(i,8);  end;

        if in_data.mater.rho~=0  rho_=in_data.mater.rho;  else  rho_=in_data.EL(i,9);     end;
        
        
        [Klc,Mlc] = D3_TETRAH (in_data.ND(node1,2), in_data.ND(node1,3), in_data.ND(node1,4),...
            in_data.ND(node2,2), in_data.ND(node2,3), in_data.ND(node2,4),...
            in_data.ND(node3,2), in_data.ND(node3,3), in_data.ND(node3,4),...
            in_data.ND(node4,2), in_data.ND(node4,3), in_data.ND(node4,4), Em, miu_1);
        
        Mlc = Mlc.*rho_;

        t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd; t4 = node4*dofPnd; 
        bg = [(t1-(dofPnd-1)):t1];  md1= [(t2-(dofPnd-1)):t2];
        md2= [(t3-(dofPnd-1)):t3];  md3= [(t4-(dofPnd-1)):t4];
        
        Kgl([bg md1 md2 md3],[bg md1 md2 md3]) = ...
            Kgl([bg md1 md2 md3],[bg md1 md2 md3]) + Klc;
        M([bg md1 md2 md3],[bg md1 md2 md3]) = ...
            M([bg md1 md2 md3],[bg md1 md2 md3]) + Mlc;
    end;
    
     if in_data.EL(i,2)==5 
         
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
        node4 = find(in_data.ND(:,1)==in_data.EL(i,6));

        if in_data.mater.E~=0  Em=in_data.mater.E;   else Em=in_data.EL(i,7);   end;
        if in_data.mater.h~=0 
            if length(in_data.mater.h)==1   hT=[1 1 1 1].*in_data.mater.h; end
            if length(in_data.mater.h)==4   hT=in_data.mater.h;            end
        else hT = in_data.EL(i,8:11); end; 
        
        if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else  miu_1 = in_data.EL(i,12); end;
        if in_data.mater.rho~=0 rho_=in_data.mater.rho;  else  rho_=in_data.EL(i,13);    end;

        if in_data.EL(i,2)==5
            [Bsys,Esys,Klc,Msys] = D2_CSQ (in_data.ND(node1,2),in_data.ND(node1,3), ...
                in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node3,2),in_data.ND(node3,3),...
                in_data.ND(node4,2),in_data.ND(node4,3), Em, hT, miu_1);
        end
        if in_data.EL(i,2)==51
        ex = in_data.PML.ex;
        ey = in_data.PML.ey;
            [Bsys,Esys,Klc,Msys] = D2_CSQ_PML (in_data.ND(node1,2),in_data.ND(node1,3), ...
                in_data.ND(node2,2),in_data.ND(node2,3),in_data.ND(node3,2),in_data.ND(node3,3),...
                in_data.ND(node4,2),in_data.ND(node4,3), Em, hT, miu_1, ex, ey);
        end
        
        M_loc   = spalloc(8,8,1);
        for z=1:8 M_loc(z,z)=1; end; 
        
        V_loc   = .5*(( in_data.ND(node3,2)-in_data.ND(node1,2) )*...
                  ( in_data.ND(node4,3)-in_data.ND(node2,3) )-...
                  ( in_data.ND(node4,2)-in_data.ND(node2,2) )*...
                  ( in_data.ND(node3,3)-in_data.ND(node1,3) ))* mean(hT);
        
        Msys = M_loc.*rho_*V_loc/4;
        if in_data.EL(i,2)==51
            Msys = Msys.*diag([ex ey ex ey ex ey ex ey]);
        end


        t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd; t4 = node4*dofPnd;
        bg = [(t1-(dofPnd-1)):t1];  md1= [(t2-(dofPnd-1)):t2];
        md2= [(t3-(dofPnd-1)):t3];  en = [(t4-(dofPnd-1)):t4];
        
        Kgl([bg md1 md2 en],[bg md1 md2 en]) = ...
            Kgl([bg md1 md2 en],[bg md1 md2 en]) + Klc;
        M([bg md1 md2 en],[bg md1 md2 en]) = ...
            M([bg md1 md2 en],[bg md1 md2 en]) + Msys;
    end;

     if in_data.EL(i,2)==9
        
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
        
        if in_data.mater.E~=0  E1=in_data.mater.E;       else E1=in_data.EL(i,6);       end;
        if in_data.mater.h~=0 h_1=in_data.mater.h;       else h_1 = in_data.EL(i,7);    end;
        if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else  miu_1 = in_data.EL(i,8); end;
        
        x1 = in_data.ND(node1,2); y1 = in_data.ND(node1,3);
        x2 = in_data.ND(node2,2); y2 = in_data.ND(node2,3);
        x3 = in_data.ND(node3,2); y3 = in_data.ND(node3,3);
        
        [Klc,ML] = D2_BCIZ (x1,y1,x2,y2,x3,y3,E1,h_1,miu_1);

        t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd;
        
        A = (0.5 * det([1 1 1; x1 x2 x3; y1 y2 y3]) );
        
        ML = ML.*in_data.mater.rhoX;
        bg = [(t1-2):t1];  md1= [(t2-2):t2]; md2= [(t3-2):t3];
        
        Kgl([bg md1 md2],[bg md1 md2]) = Kgl([bg md1 md2],[bg md1 md2]) + Klc;
        M([bg md1 md2],[bg md1 md2]) =   M([bg md1 md2],[bg md1 md2]) + ML;
    end;    
end;

sp_=sprintf('\n');
disp([sp_ '    nnz in K: ' num2str(nnz(Kgl))]);