% plot_declined_beam
 function [xord,yord,x1,y1,R,ux_new,vy_new] = frame_plot(U,i_ele,XORD,YORD,NumNP,NP,r,N)
%
% r: Deformation enlargement ratio.
% N: interpolation points in every beam element.
xord=XORD+U(1:3:3*NumNP); % new coordinates after the deformation
yord=YORD+U(2:3:3*NumNP); % 
x1=xord(NP(i_ele,1)); x2=xord(NP(i_ele,2)); % (x1,y1): the 1st node of the beam element
y1=yord(NP(i_ele,1)); y2=yord(NP(i_ele,2)); % (x2,y2): the 1st node of the beam element
theta=atan2((yord(NP(i_ele,2))-yord(NP(i_ele,1))),(xord(NP(i_ele,2))-xord(NP(i_ele,1)))); % Rotation angle of the beam element i_ele.
R=[cos(theta) -sin(theta);sin(theta) cos(theta)]; % Rotation matrix
% L0=((XORD(NP(i_ele,1))-XORD(NP(i_ele,2)))^2+(YORD(NP(i_ele,1))-YORD(NP(i_ele,2)))^2)^0.5; % Original length
L=((xord(NP(i_ele,1))-xord(NP(i_ele,2)))^2+(yord(NP(i_ele,1))-yord(NP(i_ele,2)))^2)^0.5; % Deformed length
def=[]; % Initial matrix 
i=i_ele; % the 1st beam element
phe1=U(3+(NP(i_ele,1)-1)*3); phe2=U(3+(NP(i_ele,2)-1)*3);
for n=1:N+1 % N+1 interpolation point between the 1st TWO node-points
    x=L/N*(n-1);% local coordinate    
    % get the deflection for the interpolation points in the element
    N1_x = (2*x^3-3*x^2+L^3)/L^3; % Shape function
    N2_x = (L*x^3-2*L^2*x^2+x*L^3)/L^3;
    N3_x = (-2*x^3+3*L*x^2)/L^3;
    N4_x = (L*x^3-L^2*x^2)/L^3;
    v_def = [N1_x N2_x N3_x N4_x]*[0;phe1;0;phe2]; % v_def: deflection in Local coordinates
    def=[def;x v_def];
end
ux=def(:,1)'; % x-coordinate
vy=r*def(:,2)'; % deflection in LOCAL coordinates
ux_vy_new=[]; % 
for i=1:N+1
    [ux_vy_new]=[ux_vy_new R*[ux(i);vy(i)]]; % the ux, vy in GLOBAL coordinates
end
ux_new=x1+ux_vy_new(1,:); % Transportation plus Rotation
vy_new=y1+ux_vy_new(2,:);
% figure, 
 plot(ux_new,vy_new,'r'); grid on % plot declined deformation of the beam element
hold on
plot(ux_new(1),vy_new(1),'bo'); hold on % Plot the 1st node of the deformed beam element
plot(ux_new(N+1),vy_new(N+1),'bo'); hold on % Plot the 2nd node of the deformed beam element