% -------------------------------- %
%   Milad Telecommunication Tower  %
%          (FEM Analysis)          %
%      by                          %
%        Bahman Aboulhasanzadeh    %
%        Sharif Uni. of Tech.      %
%                                  %
%          Summer - 2004           %
% -------------------------------- %

clc;
clear all;

Inputs = xlsread('D:\Documents and Settings\Bahman\Desktop\Finite Element\Full Computer Approach\inputs.xls');
for i = 1:18 
    x = Inputs(i,7)*pi/180; y= Inputs(i,8)*pi/180; z=Inputs(i,9)*pi/180; A = Inputs(i,6) ; Iy = Inputs(i,11) ; Iz = Inputs(i,12) ; E = Inputs(i,10); L =Inputs(i,5) ;
    T =[[cos(z)*cos(y),sin(z)*cos(y),-sin(y),0,0,0,0,0,0,0];[cos(z)*sin(y)*sin(x)-cos(x)*sin(z),sin(y)*sin(x)*sin(z)+cos(x)*cos(z),sin(x)*cos(y),0,0,0,0,0,0,0];[cos(x)*sin(y)*cos(z)+sin(x)*sin(z),cos(x)*sin(y)*sin(z)-sin(x)*cos(z),cos(x)*cos(y),0,0,0,0,0,0,0];[0,0,0,1,0,0,0,0,0,0];[0,0,0,0,1,0,0,0,0,0];[0,0,0,0,0,cos(z)*cos(y),sin(z)*cos(y),-sin(y),0,0];[0,0,0,0,0,cos(z)*sin(y)*sin(x)-cos(x)*sin(z),sin(y)*sin(x)*sin(z)+cos(x)*cos(z),sin(x)*cos(y),0,0];[0,0,0,0,0,cos(x)*sin(y)*cos(z)+sin(x)*sin(z),cos(x)*sin(y)*sin(z)-sin(x)*cos(z),cos(x)*cos(y),0,0];[0,0,0,0,0,0,0,0,1,0];[0,0,0,0,0,0,0,0,0,1]];
    Kl_f = [[A*E/L,0,0,0,0,-A*E/L,0,0,0,0];[0,12*E*Iz/L^3,0,0,6*E*Iz/L^2,0,-12*E*Iz/L^3,0,0,6*E*Iz/L^2];[0,0,12*E*Iy/L^3,-6*E*Iy/L^2,0,0,0,-12*E*Iy/L^3,-6*E*Iy/L^2,0];[0,0,-6*E*Iy/L^2,4*E*Iy/L,0,0,0,6*E*Iy/L^2,2*E*Iy/L,0];[0,6*E*Iz/L^2,0,0,4*E*Iz/L,0,-6*E*Iz/L^2,0,0,2*E*Iz/L];[-A*E/L,0,0,0,0,A*E/L,0,0,0,0];[0,-12*E*Iz/L^3,0,0,-6*E*Iz/L^2,0,12*E*Iz/L^3,0,0,-6*E*Iz/L^2];[0,0,-12*E*Iy/L^3,6*E*Iy/L^2,0,0,0,12*E*Iy/L^3,6*E*Iy/L^2,0];[0,0,-6*E*Iy/L^2,2*E*Iy/L,0,0,0,6*E*Iy/L^2,4*E*Iy/L,0];[0,6*E*Iz/L^2,0,0,2*E*Iz/L,0,-6*E*Iz/L^2,0,0,4*E*Iz/L]];
    K{i} = T'*Kl_f*T;
end
K_total_Big = zeros(125)
for i = 1:18
    Start_Node = Inputs(i,13)
    End_Node = Inputs(i,15)
    Index_Start = (Start_Node-1)*5
    Index_End = (End_Node-1)*5
    B = [Index_Start+1 , Index_Start+2 , Index_Start+3 , Index_Start+4 , Index_Start+5 , Index_End+1 , Index_End+2 , Index_End+3 , Index_End+4 , Index_End+5]
    K_total_Big = total_mat(K{i},B,K_total_Big)
end
 
wk1write('D:\Documents and Settings\Bahman\Desktop\Finite Element\Full Computer Approach\BigMatrix.xls',K_total_Big,1,1); 


% ----- Matrix Reduction Part
K_Cut = K_total_Big;
Dof = zeros(25,5);   % A Matrix to determine Which DOF of Nodes Constrainted (Its row & column deleted)

for i = 25:-1:1
    if Inputs(i,20)==0 || Inputs(i,20)==1
        if Inputs(i,20)==0  % Case: Node of Link
            size_K_Cut = size(K_Cut);
            K_Cut = [K_Cut(:,1:i*5-1),K_Cut(:,i*5+1:size_K_Cut(1,2))];    % Cutting t_zn Column
            K_Cut = [K_Cut(1:i*5-1,:);K_Cut(i*5+1:size_K_Cut(1,1),:)];    % Cutting t_zn Row
            Dof(i,5)=1;
            size_K_Cut = size(K_Cut);
            K_Cut = [K_Cut(:,1:i*5-2),K_Cut(:,i*5:size_K_Cut(1,2))];    % Cutting t_yn Column
            K_Cut = [K_Cut(1:i*5-2,:);K_Cut(i*5:size_K_Cut(1,1),:)];    % Cutting t_yn Row
            Dof(i,4)=1;
        else  % Case: Node of Beam
            size_K_Cut = size(K_Cut);
            K_Cut = [K_Cut(:,1:i*5-2),K_Cut(:,i*5:size_K_Cut(1,2))];    % Cutting t_yn Column
            K_Cut = [K_Cut(1:i*5-2,:);K_Cut(i*5:size_K_Cut(1,1),:)];    % Cutting t_yn Row
            Dof(i,4)=1;
        end;
        if Inputs(i,19) ~= 0   % Connected Beam to Link Nodes 
            if Inputs(i,21) == 0   % Node have Sym. Condition
                K_Cut(:,i*5-4) = K_Cut(:,i*5-4) + 0.198912 * K_Cut(:,i*5-2);     % Adding Constraint of Sym. Condition for u_n
                K_Cut(i*5-4,:) = K_Cut(i*5-4,:) + 0.198912 * K_Cut(i*5-2,:);
            end;
            K_Cut(:,Inputs(i,19)*5-2) = K_Cut(:,Inputs(i,19)*5-2) + K_Cut(:,i*5-2);   % Adding Link Beam Node Condition for w_m
            K_Cut(Inputs(i,19)*5-2,:) = K_Cut(Inputs(i,19)*5-2,:) + K_Cut(i*5-2,:);
            size_K_Cut = size(K_Cut);
            K_Cut = [K_Cut(:,1:i*5-3),K_Cut(:,i*5-1:size_K_Cut(1,2))];    % Cutting w_n Column
            K_Cut = [K_Cut(1:i*5-3,:);K_Cut(i*5-1:size_K_Cut(1,1),:)];    % Cutting w_n Row
            Dof(i,3)=1;
            K_Cut(:,Inputs(i,19)*5-3) = K_Cut(:,Inputs(i,19)*5-3) + K_Cut(:,i*5-3);   % Adding Link Beam Node Condition for v_m
            K_Cut(Inputs(i,19)*5-3,:) = K_Cut(Inputs(i,19)*5-3,:) + K_Cut(i*5-3,:);
            size_K_Cut = size(K_Cut);
            K_Cut = [K_Cut(:,1:i*5-4),K_Cut(:,i*5-2:size_K_Cut(1,2))];    % Cutting v_n Column
            K_Cut = [K_Cut(1:i*5-4,:);K_Cut(i*5-2:size_K_Cut(1,1),:)];    % Cutting v_n Row
            Dof(i,2)=1;
            K_Cut(:,Inputs(i,19)*5-4) = K_Cut(:,Inputs(i,19)*5-4) + K_Cut(:,i*5-4);   % Adding Link Beam Node Condition for u_m
            K_Cut(Inputs(i,19)*5-4,:) = K_Cut(Inputs(i,19)*5-4,:) + K_Cut(i*5-4,:);
            size_K_Cut = size(K_Cut);
            K_Cut = [K_Cut(:,1:i*5-5),K_Cut(:,i*5-3:size_K_Cut(1,2))];    % Cutting u_n Column
            K_Cut = [K_Cut(1:i*5-5,:);K_Cut(i*5-3:size_K_Cut(1,1),:)];    % Cutting u_n Row
            Dof(i,1)=1;
        else
            if Inputs(i,21) == 0   % Node have Sym. Condition
                K_Cut(:,i*5-4) = K_Cut(:,i*5-4) + 0.198912 * K_Cut(:,i*5-2);     % Adding Constraint of Sym. Condition for u_n
                K_Cut(i*5-4,:) = K_Cut(i*5-4,:) + 0.198912 * K_Cut(i*5-2,:);
                size_K_Cut = size(K_Cut);
                K_Cut = [K_Cut(:,1:i*5-3),K_Cut(:,i*5-1:size_K_Cut(1,2))];    % Cutting w_n Column
                K_Cut = [K_Cut(1:i*5-3,:);K_Cut(i*5-1:size_K_Cut(1,1),:)];    % Cutting w_n Row
                Dof(i,3)=1;
            elseif Inputs(i,21) == 1
                size_K_Cut = size(K_Cut);
                K_Cut = [K_Cut(:,1:i*5-3),K_Cut(:,i*5-1:size_K_Cut(1,2))];    % Cutting w_n Column
                K_Cut = [K_Cut(1:i*5-3,:);K_Cut(i*5-1:size_K_Cut(1,1),:)];    % Cutting w_n Row
                Dof(i,3)=1;
            end;
        end;
    elseif Inputs(i,20)==2
        size_K_Cut = size(K_Cut);
        K_Cut = [K_Cut(:,1:i*5-1),K_Cut(:,i*5+1:size_K_Cut(1,2))];    % Cutting t_zn Column
        K_Cut = [K_Cut(1:i*5-1,:);K_Cut(i*5+1:size_K_Cut(1,1),:)];    % Cutting t_zn Row
        Dof(i,5)=1;
        size_K_Cut = size(K_Cut);
        K_Cut = [K_Cut(:,1:i*5-2),K_Cut(:,i*5:size_K_Cut(1,2))];    % Cutting t_yn Column
        K_Cut = [K_Cut(1:i*5-2,:);K_Cut(i*5:size_K_Cut(1,1),:)];    % Cutting t_yn Row
        Dof(i,4)=1;
        size_K_Cut = size(K_Cut);
        K_Cut = [K_Cut(:,1:i*5-3),K_Cut(:,i*5-1:size_K_Cut(1,2))];    % Cutting w_n Column
        K_Cut = [K_Cut(1:i*5-3,:);K_Cut(i*5-1:size_K_Cut(1,1),:)];    % Cutting w_n Row
        Dof(i,3)=1;
        size_K_Cut = size(K_Cut);
        K_Cut = [K_Cut(:,1:i*5-4),K_Cut(:,i*5-2:size_K_Cut(1,2))];    % Cutting v_n Column
        K_Cut = [K_Cut(1:i*5-4,:);K_Cut(i*5-2:size_K_Cut(1,1),:)];    % Cutting v_n Row
        Dof(i,2)=1;    
        size_K_Cut = size(K_Cut);
        K_Cut = [K_Cut(:,1:i*5-5),K_Cut(:,i*5-3:size_K_Cut(1,2))];    % Cutting u_n Column
        K_Cut = [K_Cut(1:i*5-5,:);K_Cut(i*5-3:size_K_Cut(1,1),:)];    % Cutting u_n Row
        Dof(i,1)=1;   
    end;
end;
%  template = {'u','v','w','ty','tz'};               
horizontal_title = [];
vertical_title = [];
for i= 1:25
     for j= 1:5
         if Dof(i,j) == 0 
             horizontal_title = [horizontal_title,str2num(strcat(num2str(j),'00',num2str(i)))];
             vertical_title = [vertical_title;str2num(strcat(num2str(j),'00',num2str(i)))];
         end;
     end;
end;
 
wk1write('D:\Documents and Settings\Bahman\Desktop\Finite Element\Full Computer Approach\Totalmatrix_Cut.xls',K_Cut,1,1); 
wk1write('D:\Documents and Settings\Bahman\Desktop\Finite Element\Full Computer Approach\horizontal_title.xls',horizontal_title,1,1);         
wk1write('D:\Documents and Settings\Bahman\Desktop\Finite Element\Full Computer Approach\vertical_title.xls',vertical_title,1,1);