Code covered by the BSD License  

Highlights from
FEM toolbox for solid mechanics

image thumbnail
from FEM toolbox for solid mechanics by Anton Zaicenco
The finite element toolbox for solid mechanics with GUI.

plot_2dbeam_stat_diagram (in_data,obj,resp,type,fign)
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');

Contact us at files@mathworks.com