k_assembly.m

FEM tools for caculation of nonlinear problems

function [Kgl, M] = K_assembly (in_data, type_analysis, part_N_r, part_N_c)

    [dofPnd,EL_TYPE] = type_of_elem (in_data);
    
    dof = size(in_data.ND)*dofPnd;
    nn_z = floor((1+part_N_r(2)-part_N_r(1))*(1+part_N_c(2)-part_N_c(1))*0.05);

   Kgl = spalloc(1+part_N_r(2)-part_N_r(1), 1+part_N_c(2)-part_N_c(1),nn_z);

    if  isfield(in_data,'T')
        M = spalloc(1+part_N_r(2)-part_N_r(1), 1,nn_z);
    else
        M   = spalloc(1+part_N_r(2)-part_N_r(1), 1+part_N_c(2)-part_N_c(1),nn_z);
    end

indj =1;


sxEL = size(in_data.EL,1);

for i=1:sxEL

    if i>(indj*250) 
        disp(['     K assembly: ' num2str(indj*250) ' elements ']); 
        indj=indj+1; 
    end;
    if i==sxEL 
        disp(['     K assembly: ' num2str(i) ' elements - complete']);  
    end;
    
     if in_data.EL(i,2)==44
         
         dofPnd = 1;
         
         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));
         a=in_data.T.ND;
         in_data.T.Tmatrix = zeros(size(in_data.EL,1),3);
         in_data.T.hmatrix = zeros(size(in_data.EL,1),3);
         for i=1:size(in_data.EL,1)
             b1 = in_data.EL(i,3:4);
             b2 = in_data.EL(i,4:5);
             b3 = [in_data.EL(i,5)  in_data.EL(i,3)];
             ind1=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(1) ones(1,length(a)-1)'*b1(2)]')));
             ind1_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(2) ones(1,length(a)-1)'*b1(1)]')));
             ind2=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(1) ones(1,length(a)-1)'*b2(2)]')));
             ind2_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(2) ones(1,length(a)-1)'*b2(1)]')));
             ind3=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(1) ones(1,length(a)-1)'*b3(2)]')));
             ind3_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(2) ones(1,length(a)-1)'*b3(1)]')));
             if [ind1 ind1_]
                 in_data.T.Tmatrix(i,1) = in_data.T.T([ind1 ind1_]); 
                 in_data.T.hmatrix(i,1)= in_data.T.h([ind1 ind1_]);
             end
             if [ind2 ind2_]
                 in_data.T.Tmatrix(i,2) = in_data.T.T([ind2 ind2_]); 
                 in_data.T.hmatrix(i,2)= in_data.T.h([ind2 ind2_]);
             end
             if [ind3 ind3_]
                 in_data.T.Tmatrix(i,3) = in_data.T.T([ind3 ind3_]); 
                 in_data.T.hmatrix(i,3)= in_data.T.h([ind3 ind3_]);
             end
         end

         a=in_data.flux.ND;
         in_data.flux.qmatrix = zeros(size(in_data.EL,1),3);
         for i=1:size(in_data.EL,1)
             b1 = in_data.EL(i,3:4);
             b2 = in_data.EL(i,4:5);
             b3 = [in_data.EL(i,5)  in_data.EL(i,3)];
             ind1=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(1) ones(1,length(a)-1)'*b1(2)]')));
             ind1_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b1(2) ones(1,length(a)-1)'*b1(1)]')));
             ind2=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(1) ones(1,length(a)-1)'*b2(2)]')));
             ind2_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b2(2) ones(1,length(a)-1)'*b2(1)]')));
             ind3=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(1) ones(1,length(a)-1)'*b3(2)]')));
             ind3_=find(~sum(abs([a([1:end-1]); a(2:end)] - [ones(1,length(a)-1)'*b3(2) ones(1,length(a)-1)'*b3(1)]')));
             if [ind1 ind1_]
                 in_data.flux.qmatrix(i,1) = in_data.flux.q([ind1 ind1_]); 
             end
             if [ind2 ind2_]
                 in_data.flux.qmatrix(i,2) = in_data.flux.q([ind2 ind2_]); 
             end
             if [ind3 ind3_]
                 in_data.flux.qmatrix(i,3) = in_data.flux.q([ind3 ind3_]); 
             end
         end
         
         kx = in_data.EL(i,6);
         ky = in_data.EL(i,7);
         Q  = in_data.EL(i,8);
         h  = in_data.T.hmatrix(i,:);
         Tinf = in_data.T.Tmatrix(i,:);
         q  = in_data.flux.qmatrix(i,:);
         
         
         [Klc,rlc] = D2_CST_temp (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), ...
             kx, ky, h, Q, Tinf, q);
         
         t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd;
         
         bg = [(t1-(dofPnd-1)):t1]; md = [(t2-(dofPnd-1)):t2];
         en = [(t3-(dofPnd-1)):t3];
         
         Kgl([bg md en],[bg md en]) = Kgl([bg md en],[bg md en]) + Klc;
         M([bg md en]) = M([bg md en]) + rlc;
         
         
     end

    if in_data.EL(i,2)==0 | in_data.EL(i,2)==1 | in_data.EL(i,2)==2 | in_data.EL(i,2)==31

        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));

        if in_data.mater.E~=0    E1=in_data.mater.E;      else E1=in_data.EL(i,5);     end;
        if in_data.mater.A~=0    A_1=in_data.mater.A;     else A_1 = in_data.EL(i,6);  end;
        if in_data.mater.I~=0    I_1=in_data.mater.I;     else I_1 = in_data.EL(i,7);  end;
        if in_data.mater.rho~=0  rho_=in_data.mater.rho;  else rho_= in_data.EL(i,8);  end;

        if in_data.EL(i,2)==0
            [M_loc,Klc] = FF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),...
                in_data.ND(node2,2), in_data.ND(node2,3),E1,A_1,I_1,rho_);
        end;
        if in_data.EL(i,2)==1
            [M_loc,Klc] = FP_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),...
                in_data.ND(node2,2), in_data.ND(node2,3),E1,A_1,I_1,rho_);
        end;
        if in_data.EL(i,2)==2
            [M_loc,Klc] = PF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),...
                in_data.ND(node2,2), in_data.ND(node2,3),E1,A_1,I_1,rho_);
        end;
        if in_data.EL(i,2)==31
            [M_loc,Klc] = D3_TRUSS (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, A_1, rho_);
        end;
        
        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)==4
        
        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;
        if in_data.mater.rho~=0  rho_=in_data.mater.rho;  else rho_= in_data.EL(i,9);    end;
        

        
        [Bsys,Esys,Klc,Msys] = D2_CST (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), E1,h_1,miu_1);

        Msys = Msys.*rho_;
        
        t1 = node1*dofPnd; t2 = node2*dofPnd; t3 = node3*dofPnd;

        bg = [(t1-(dofPnd-1)):t1]; md = [(t2-(dofPnd-1)):t2];
        en = [(t3-(dofPnd-1)):t3];
        
        Kgl([bg md en],[bg md en]) = Kgl([bg md en],[bg md en]) + Klc;
        M([bg md en],[bg md en]) = M([bg md en],[bg md en]) + Msys;
    end;
 
     if in_data.EL(i,2)==3
        
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));