function U = fem(E,P)
%% Deflection of a beam using a 4 node isoparametric element
% Coded by(Matthias Schubert and Donat Perler, 2007, PhD Seminar 'The Finite Element Method and the Analysis of Systems
% with Uncertain Properties')

%% Input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Geometrical input 
l=1000;          % Length   [mm]
h=60;            % Height   [mm]
b_el=40;         % Width    [mm]

Num_el_x=30;        % Number of elements in x-direction
Num_el_y=10;         % Number of elements in y-direction

Upper_nodes=1:(Num_el_y+1):(Num_el_x+1)*(Num_el_y+1)-Num_el_y;
Lower_nodes=(Num_el_y+1):(Num_el_y+1):(Num_el_x+1)*(Num_el_y+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%Input forces

% Force at node number
F_Nodes_x=[];%[1+(Num_el_y+1)*Num_el_x];                                           
F_Nodes_y=[1+(Num_el_y+1)*Num_el_x];                       

% Force
Force_x=[];          %[N]
Force_y=[P];         %[N]

% Boundary conditions at nodes
B_Node_x=[1 2 3 4 5 6 7 8 9];       %Node Number
B_Node_y=[1 2 3 4 5 6 7 8 9];       %Node Number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Material Properties
nu=0.3;     %Poisson ratio
%E=210000;   %Stiffness [N/mm^2]

% Aplification for plot
% amp=100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Main Programm
%%
% Number of Elements
    N_El=Num_el_x*Num_el_y;
%%
% Number of nodes
    N_Nodes=(Num_el_x+1)*(Num_el_y+1);
%%
% Generating the Load Vector
    Load_glob=zeros(N_Nodes*2,1);
    for i=1:size(F_Nodes_y,2)
        Load_glob(F_Nodes_y(i)*2,:)=Force_y(i);
    end
    for i=1:size(F_Nodes_x,2)
        Load_glob(F_Nodes_x(i)*2-1,:)=Force_x(i);
    end
%%
% Generating boundary vector
    for i=1:size(B_Node_x,2)
        Boundary(i,1)=B_Node_x(i)*2-1;
    end

    for j=1:size(B_Node_y,2)
        Boundary(j+i,1)=B_Node_y(j)*2;
    end
    Boundary=sort(Boundary);
%%
% Generating nodal coordinates
    X_coordinates=zeros(Num_el_y+1,Num_el_x+1);
    Y_coordinates=zeros(Num_el_y+1,Num_el_x+1);
    for x_coor=1:(Num_el_y+1)
        X_coordinates(x_coor,:)=linspace(0,l,(Num_el_x+1));
    end
    for y_coor=1:(Num_el_x+1)
        Y_coordinates(:,y_coor)=linspace(0,h,(Num_el_y+1));
    end

    trans1=(1:1:Num_el_y+1)';
    trans2=(size(trans1,1)+1)*ones(Num_el_y+1,1);
    trans=trans2-trans1;
    clear trans1 trans2;
    Y_coordinates=Y_coordinates(trans,:);
%%
% Assembling elements
    Index=[2 4 3 1];
    ele_counter=0;
    Node_tmp_1=1;
    Node_tmp_2=Num_el_y+2;

    for m=1:Num_el_x

        for n=1:Num_el_y

             ele_counter=ele_counter+1;
             eval(sprintf('Element_%i=[Node_tmp_1 Node_tmp_1+1 Node_tmp_2 Node_tmp_2+1];',ele_counter));   
             eval(sprintf('Element_%i=Element_%i(Index);',ele_counter,ele_counter));
             Node_tmp_1=Node_tmp_1+1;
             Node_tmp_2=Node_tmp_2+1;

        end    
         Node_tmp_1=Node_tmp_1+1;
         Node_tmp_2=Node_tmp_2+1;

    end
%%
% Local element stiffness matrix (analytically derived)
    a_el=l/Num_el_x;
    h_el=h/Num_el_y;
    gamma=a_el/h_el;
    Psi1=(1+nu)*gamma;
    Psi2=(1-3*nu)*gamma;
    Psi3=2+(1-nu)*gamma^2;
    Psi4=2*gamma^2+(1-nu);
    Psi5=(1-nu)*gamma^2-4;
    Psi6=(1-nu)*gamma^2-1;
    Psi7=4*gamma^2-(1-nu);
    Psi8=gamma^2-(1-nu);


    k_4node =(E*b_el)/(24*gamma*(1-nu^2))*[4*Psi3    3*Psi1      2*Psi5       -3*Psi2        -2*Psi3     -3*Psi1     -4*Psi6     3*Psi2      ;...
                                          3*Psi1    4*Psi4      3*Psi2       4*Psi8         -3*Psi1     -2*Psi4     -3*Psi2     -2*Psi7     ;...
                                          2*Psi5    3*Psi2      4*Psi3       -3*Psi1        -4*Psi6     -3*Psi2     -2*Psi3     3*Psi1      ;...
                                          -3*Psi2   4*Psi8      -3*Psi1      4*Psi4         3*Psi2      -2*Psi7     3*Psi1      -2*Psi4     ;...
                                          -2*Psi3   -3*Psi1     -4*Psi6      3*Psi2         4*Psi3      3*Psi1      2*Psi5      -3*Psi2     ;...
                                          -3*Psi1   -2*Psi4     -3*Psi2      -2*Psi7        3*Psi1      4*Psi4      3*Psi2      4*Psi8      ;...
                                          -4*Psi6   -3*Psi2     -2*Psi3      3*Psi1         2*Psi5      3*Psi2      4*Psi3      -3*Psi1     ;...
                                          3*Psi2    -2*Psi7     3*Psi1       -2*Psi4        -3*Psi2      4*Psi8      -3*Psi1     4*Psi4      ];

%%                                     
% Assemble global stiffness matrix
     K_global=sparse(N_Nodes*2,N_Nodes*2);
     for i=1:N_El
         for j=1:4 
             for k=1:4      
                 eval(sprintf('K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2-1 )= K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2-1 )  + k_4node(j*2-1,k*2-1);',i,i,i,i));     
                 eval(sprintf('K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2   )= K_global( Element_%i(1,j)*2-1, Element_%i(1,k)*2 )    + k_4node(j*2-1,k*2);',i,i,i,i));       
                 eval(sprintf('K_global( Element_%i(1,j)*2,   Element_%i(1,k)*2-1 )= K_global( Element_%i(1,j)*2, Element_%i(1,k)*2-1 )    + k_4node(j*2,k*2-1);',i,i,i,i));        
                 eval(sprintf('K_global( Element_%i(1,j)*2,   Element_%i(1,k)*2   )= K_global( Element_%i(1,j)*2, Element_%i(1,k)*2 )      + k_4node(j*2,k*2);',i,i,i,i));
             end
         end
     end
     
%% 
% Result
        K_global_tmp=K_global;
        K_global(:,Boundary)=[];
        K_global(Boundary,:)=[];
        Load_glob(Boundary,:)=[];
        UG_result=K_global\Load_glob;
%%
% Creating deformation marix
       U=(zeros(N_Nodes*2,1));
       Ug_init=(1:1:N_Nodes*2)';
       Ug_init(Boundary,:)=[];
       U(Ug_init,:)=UG_result;
 
%       U_x_tmp=(U(1:2:end).*amp)';
%       U_y_tmp=(U(2:2:end).*amp)';
% 
%       U_x=zeros(Num_el_y+1,Num_el_x+1);
%       U_y=zeros(Num_el_y+1,Num_el_x+1);
% 
%       for t=1:(Num_el_y+1)
%           eval(sprintf('U_x(%i,:)=U_x_tmp(%i:Num_el_y+1:end);',t,t));
%           eval(sprintf('U_y(%i,:)=U_y_tmp(%i:Num_el_y+1:end);',t,t));
%       end