MEFtave/Sol_analytique/poutre_compar.m

function  poutre_compar(U)
% comparaison de la solution numerique et analytique pour 
% une poutre uniformement chargee sur deux appuis
%  H.Oudin 
global nddln nelt nnode
global Coord Connec Nprop Prop

fa1=@(x) Prop(1,2)*x.*(x-1)/2 ;  % solution analytique
fa2=@(x) -Prop(1,2)*(2*x-1)/2;
taille = get(0,'ScreenSize'); 
figure('Name','comparaison avec la solution analytique',...
      'Position',[taille(3)/2.02 taille(4)/2.6 taille(3)/2 taille(4)/2]) 

subplot(2,1,1),title('moment de flexion Mf'), hold on,     
fplot(fa1,[0 1],'r')   %----- solution analytique
for iel=1:nelt          %----- boucle sur les elements solution numerique
  loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]]; end
	Ue=U(loce);
  EI=Prop(Nprop(iel),1);
  X = Coord(Connec(iel,:)); x2 = X(2); x1=X(1); L = abs(x2-x1); 
  Mfe= (EI/L^2)*[-6*Ue(1)-4*L*Ue(2)+6*Ue(3)-2*L*Ue(4),... 
                  6*Ue(1)+2*L*Ue(2)-6*Ue(3)+4*L*Ue(4)];
  x=x1:(x2-x1)/5:x2; y=Mfe(1)+(Mfe(2)-Mfe(1))*(x-x1)/L;
  plot(x,y,'b')
end
grid
subplot(2,1,2)
hold on,       
fplot(fa2,[0 1],'r')
title('effort tranchant'),
for iel=1:nelt          %----- boucle sur les elements solution numerique
    loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
	Ue=U(loce);
    EI=Prop(Nprop(iel),1);
    X = Coord(Connec(iel,:)); x2 = X(2); x1=X(1); L = abs(x2-x1); 
    Te = -(EI/L^3)*(12*Ue(1)+6*L*Ue(2)-12*Ue(3)+6*L*Ue(4));
    x=[x1:(x2-x1)/5:x2];
    y=Te*ones(1,length(x));
    plot(x,y,'b'),
end
grid
return