Plotting deformation of frame structure (FEM)

This is a MATLAB code to plot a deformation of planar frame structure.

Since the formulation involves transformation of local elements to a global system,

plotting using a simple w(x) versus x is not possible..

We need to transform back the DOFs obtained in the global system to the local system,

and plot in a local system and then simply ‘move & rotate’ it to suit the global position..

For example:

—————————————————————-

clear
clc

% INPUT 
Points = [0 0; 5 0; 8 -6];        %Nodes
Cn = [1 2; 2 3];                  %connectivity matrix
ne = length(Cn);                  %no.of elements
U = [0 0 0 0 0 -0.8 1.5 0 0.7]';  %deformation in global system

% PLOTTING
clf
hold on

% Simple plot for intial & deformed shapes
Uplot = U; Uplot(3:3:end) = [];
NewPoints = Points + reshape(Uplot,size(Points'))';
for i = 1:ne
    plot(Points(Cn(i,:),1), ...
        Points(Cn(i,:),2),'bo--')     %initial
    plot(NewPoints(Cn(i,:),1),...
        NewPoints(Cn(i,:),2),'go-.')  %deformed
end

for i = 1:ne
    
    % Length of element 
    X = Points(Cn(i,2),1)-Points(Cn(i,1),1);
    Y = Points(Cn(i,2),2)-Points(Cn(i,1),2);
    L = sqrt(X^2+Y^2);
    
    % Angle & transformation matrix
    B = atand(Y/X);
    c = cosd(B);
    s = sind(B);
    T = [c s 0 0 0 0; -s c 0 0 0 0;
        0 0 1 0 0 0; 0 0 0 c s 0;
        0 0 0 -s c 0; 0 0 0 0 0 1];
    
    % Index number in the global system
    index = [3*Cn(i,1)-2:3*Cn(i,1)...
        3*Cn(i,2)-2:3*Cn(i,2)];
    
    u = T*U(index);  %local dofs
    x = 0:L/20:L;    %original local axis
    
    % local horizontal displacements function, u(x)
    N1bar = 1-x/L;
    N2bar = x/L;
    xx = x + N1bar*u(1) + N2bar*u(4);
    
    % local vertical displacements function, w(x)
    N1beam = -3*x.^2/L^2 + 1 + 2*x.^3/L^3;
    N2beam = -2*x.^2/L + x + x.^3/L^2;
    N3beam = 3*x.^2/L^2 - 2*x.^3/L^3;
    N4beam = -x.^2/L + x.^3/L^2;
    yy = N1beam*u(2) + N2beam*u(3) + N3beam*u(5) + N4beam*u(6);

    % adjustment to the global axis (move the origin + rotate)
    xr = Points(Cn(i,1),1) + xx*cosd(B) - yy*sind(B);
    yr = Points(Cn(i,1),2) + xx*sind(B) + yy*cosd(B);

    %plot(xx,yy,'c-') %plot at local axis (to check)
    plot(xr,yr,'k-','LineWidth',3)  %plot at global axis  
    
end
hold off
axis equal