%% One dimensional Finite Element Method
%% June 12, 2008 
%% MATH493 (Summer 2008)
%% Axial deformation of a bar

function femaxialbar

    %% Definitions 
    %%(NOTE: THIS MUST BE CHANGED TO ACCOMODATE FUNCTIONS)
    A = 1.0;
    E = 1.0;
    L = 1.0;
    R = 1.0;
    q = 1.0;
    
    
    %% Element description
    e = 3;      % number of elements
    h = L/e;    % element size
    n = e+1;    % number of nodes
    
    
    %% Node description
    for i=1:n
        node(i)=(i-1)*h;
    end
    
    
    %% Connectivity (Pointer) Matrix
    for i=1:e
        elem(i,:)=[i i+1];
    end
    
    
    %% Initializing the Stiffness Matrix and the Load Vector
    K = zeros(n);
    f = zeros(n,1);
    
    
    %% Creating the K and f (with assembly using connectivity matrix)    
    for i=1:e
        node1 = elem(i,1);      
        node2 = elem(i,2);  
        Ke = elementStiffness(A, E, h);
        fe = elementLoad(node(node1),node(node2),q,h);
        K(node1:node2,node1:node2) = K(node1:node2,node1:node2) + Ke;
        f(node1:node2)=f(node1:node2) + fe;
    end
    
    
    %% Concentrated load at right end point
    f(n) = f(n) + R;
    
    
    %% Condensation and System solve
    %% (NOTE: THIS MUST BE CHANGED AS PER THE BOUNDARY CONDITION PRESCRIBED)
    Kred = K(2:n,2:n);
    fred = f(2:n)  
    c = Kred\fred;
    
    
    %% Solution (coefficients)
    csol = [0 c'];

    
    %% Plot of the finite element solution
    %% (NOTE: THIS MUST ALSO SHOW THE EXACT SOLUTION)
    plot(node, csol, 'ro--', 'LineWidth', 3);
    
    
    %% Element Stiffness Matrix
    %% (NOTE: THIS MUST BE MODIFIED TO ACCOUNT FOR VARYING A(x),E(x)
    %%  AND COMPUTATIONS OVER THE REFERENCE ELEMENT WITH GAUSS QUADRATURE)
    function [Ke] = elementStiffness(A, E, h)
        Ke = (A*E/h)*[1 -1;-1 1];

        
    %% Element Load Vector 
    %% (NOTE: THIS MUST BE MODIFIED TO ACCOUNT FOR VARYING q(x)
    %%  AND COMPUTATIONS OVER THE REFERENCE ELEMENT WITH GAUSS QUADRATURE)
    function [fe] = elementLoad(node1,node2,q,h)
        x1 = node1;
        x2 = node2;
        fe1 = q*x2/(2*h)*(x2^2-x1^2) - q/(3*h)*(x2^3-x1^3);
        fe2 = -q*x2/(2*h)*(x2^2-x1^2) + q/(3*h)*(x2^3-x1^3);
        fe = [fe1;fe2];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        
     