plot_2dbeam_stat_diagram.m

FEM tools for caculation of nonlinear problems

function [NVM] = plot_2dbeam_stat_diagram (in_data,obj,resp,type,fign)

% PLOT STATIC DISPLACEMENTS OF 2D BEAMS

dof = size(in_data.ND,1)*3; % total dof
maxX = max(in_data.ND(:,2)); minX = min(in_data.ND(:,2));
maxY = max(in_data.ND(:,3)); minY = min(in_data.ND(:,3));
labx = (maxX / 6); laby = (maxY / 6);
labx = min([labx laby]); laby = labx;

figure(fign);
hold on;

xL = [0:0.05:1]';		
N  = length(xL);
oo = ones(N,1);
maxM=-1e50; maxN=-1e50; maxV=-1e50;

if type == 'u'  hold off; return;          end
if type == 'd'  title ('deformation');     end

if type == 'n'  h=title ({'axial forces: fraction of {\bf N}_{max}'});    end
if type == 'v'  h=title ({'shear forces: fraction of {\bf V}_{max}'});    end
if type == 'm'  h=title ({'bending moments: fraction of {\bf M}_{max}'}); end


for b=1:size(in_data.EL,1)
    node1 = find(in_data.ND(:,1)==in_data.EL(b,3));
    node2 = find(in_data.ND(:,1)==in_data.EL(b,4));
    x1 = in_data.ND(node1,2);  y1 = in_data.ND(node1,3);	
    x2 = in_data.ND(node2,2);  y2 = in_data.ND(node2,3);	
    L = sqrt( (x2-x1)^2 + (y2-y1)^2 ) ;	    
    
    if in_data.mater.E~=0  E1=in_data.mater.E;  else E1=in_data.EL(b,5);     end;
    if in_data.mater.A~=0  A_1=in_data.mater.A; else A_1 = in_data.EL(b,6);  end;
    if in_data.mater.I~=0  I_1=in_data.mater.I; else I_1 = in_data.EL(b,7);  end;
    if in_data.mater.rho~=0  rho_=in_data.mater.rho;  else rho_= in_data.EL(b,8);    end;

    if in_data.EL(b,2)==0
        [Mlc,Klc] = FF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),...
            in_data.ND(node2,3),E1,A_1,I_1,rho_);
    end;
    if in_data.EL(b,2)==1
        [Mlc,Klc] = FP_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),...
            in_data.ND(node2,3),E1,A_1,I_1);
    end;
    if in_data.EL(b,2)==2
        [Mlc,Klc] = PF_BEAM (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),...
            in_data.ND(node2,3),E1,A_1,I_1);
    end;
    g = [ 3*node1-2 ; 3*node1-1 ; 3*node1 ; 3*node2-2 ; 3*node2-1 ; 3*node2 ];
    F(:,b) = Klc * resp.static.D(g)';
    Ne = [ -F(1,b) ; F(4,b) ];	
    Ve = [ -F(2,b) ; F(5,b) ];
    Me = [ -F(3,b) ; F(6,b) ];
    
    nf = Ne(1) + [0:N-1]'*(Ne(2)-Ne(1))/(N-1);	
    vf = Ve(1) + [0:N-1]'*(Ve(2)-Ve(1))/(N-1);	
    
    mf = Me(1) + [0:N-1]'*(Me(2)-Me(1))/(N-1) +  ...
	0.5*in_data.q(b)*[0:N-1]'.*[-N+1:0]' * L^2/N^2 ;	
    
    Np(:,b) = nf;
    Vp(:,b) = vf;
    Mp(:,b) = mf;
    
    NVM (b,:) = [b Ne' Ve' Me'];
    
    if maxN<max(max(abs(nf))) maxN=max(max(abs(nf))); aNm=b; end 
    if maxV<max(max(abs(vf))) maxV=max(max(abs(vf))); sVm=b; end 
    if maxM<max(max(abs(mf))) maxM=max(max(abs(mf))); bMm=b; end 
end;

max_V=max(max(Vp)); max_N=max(max(Np)); max_M=max(max(Mp));
min_V=min(min(Vp)); min_N=min(min(Np)); min_M=min(min(Mp));
for i=1:size(Vp,1)
    Vpn(i,:) = Vp(i,:)./maxV;
    Npn(i,:) = Np(i,:)./maxN;
    Mpn(i,:) = Mp(i,:)./maxM;
end;

for b=1:size(in_data.EL,1)	
    node1 = find(in_data.ND(:,1)==in_data.EL(b,3));
    node2 = find(in_data.ND(:,1)==in_data.EL(b,4));
    x1 = in_data.ND(node1,2);  y1 = in_data.ND(node1,3);
    x2 = in_data.ND(node2,2);  y2 = in_data.ND(node2,3);
    L = sqrt( (x2-x1)^2 + (y2-y1)^2 ) ;	 
    c = (x2-x1) / L;	
    s = (y2-y1) / L;	
    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 ];

    v = [0;0;0;0;0;0];
    u = T*v;		
    ld = L + u(4) - u(1);	
    
    x = u(1) + xL*ld;
    xd =  [ x  x.*0 ] * [ c s ; -s  c ]; 

    if type == 'd'
     plot(x1+xd(:,1),y1+xd(:,2),'r-');	
    end
    if type == 'n'
     xn =  [ x Npn(:,b)*labx ] * [ c s ; -s  c ];
     plot(x1+xn(:,1),y1+xn(:,2),'r-');   
     if b==aNm
         gf=text(x1+xn(3,1),y1+xn(3,2),num2str(.01*ceil(Npn(1,b)*100)));       set_font(gf);
         gf=text(x1+xn(18,1),y1+xn(18,2),num2str(.01*ceil(Npn(end,b)*100) ));  set_font(gf);
     end;
     plot(x1+xd(:,1),y1+xd(:,2),'k-');	
     set(fign,'name',['Max =  '  num2str(max_N) ...
             '   Min = '  num2str(min_N)],'NumberTitle','off');
    end
    if type == 'v'
     xv =  [ x Vpn(:,b)*labx ] * [ c s ; -s  c ];
     plot(x1+xv(:,1),y1+xv(:,2),'r-');   
     if b==sVm
         gf=text(x1+xv(3,1),y1+xv(3,2),num2str(.01*ceil(Vpn(1,b)*100)));       set_font(gf);
         gf=text(x1+xv(18,1),y1+xv(18,2),num2str( .01*ceil(Vpn(end,b)*100) )); set_font(gf);
     end
     plot(x1+xd(:,1),y1+xd(:,2),'k-');	
     set(fign,'name',['Max =  '  num2str(max_V) ...
             '   Min = '  num2str(min_V)],'NumberTitle','off');
    end
    if type == 'm'
     xm =  [ x Mpn(:,b)*1*labx ] * [ c s ; -s  c ];
     plot(x1+xm(:,1),y1+xm(:,2),'r-');   
     if b==bMm
         gf=text(x1+xm(1,1),y1+xm(1,2),num2str(.01*ceil(Mpn(1,b)*100)));       set_font(gf);
         gf=text(x1+xm(10,1),y1+xm(10,2),num2str(.01*ceil(Mpn(10,b)*100)));    set_font(gf);
         gf=text(x1+xm(21,1),y1+xm(21,2),num2str( .01*ceil(Mpn(end,b)*100) )); set_font(gf);
     end
     plot(x1+xd(:,1),y1+xd(:,2),'k-');	 
     set(fign,'name',['Max =  '  num2str(max_M) ...
             '   Min = '  num2str(min_M)],'NumberTitle','off');
    end
end

axis off; axis equal;



% =======================================================
function set_font (gf)

set(gf,'FontSize',7,'Color','b');