Original MEFtave from Herve Oudin

Data/NUM_exo1.m

@@ -0,0 +1,86 @@
+close all ; clc;
+%  Exercice 1 du chapitre de cours NUM - elements T3
+%  Maillage en 2 element d'une plaque rectangulaire
+%
+%  Fonctions utilisees
+%   statique        : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec numero noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   T3_ep_stress    : calcul de la contrainte dans un element T3
+%   plot_sig        : visualisation de l'etat de contrainte
+%
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('Maillage en deux T3 de plaque de l''exercice 1 du chapitre NUM');
+disp('===================');
+% definition du maillage
+L = 2 ; h = 1;
+Coord=[ 0  0; L  0; L  h;0  h ];   % definition des coordonnees des noeuds
+[nnod,ndim]=size(Coord);    
+nddln=2;  nddlt=nddln*nnod;  
+Connec=[ 1 2 4; ...            % Tableau de connectivite i , j
+         2 3 4 ];
+[nelt,nnode]=size(Connec);
+Typel = 'T3_ep';               % definition du type des elements
+for i=1:nelt   Typel = char('T3_ep',Typel); end
+
+% definition des caracteristiques mecaniques
+Nprop=[1;1;1;1];           
+Prop=[ 210000 0.3 2 0 0 ];     % valeurs  (E nu e fx fy)
+
+% CL en deplacement                                     
+CL=[ 1 , 1 , 1; ...            
+     2 , 0 , 1; ...
+     4 , 1 , 1;];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);                        
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+% charges nodales: numero du noeud , Fx,Fy
+Charg=[ 3  5  0 ];              
+F=zeros(nddlt,1);	   
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+ end
+[Fx,Fy,Fz] = feval('resultante',F);  %----- resultante des charges nodales
+
+plotstr  % trace du maillage 
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statique;   % ----- resolution du probleme (pivot = 1)
+%-----	post-traitement
+plotdef(U)
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)                                      
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+    loce=[]; 
+    for i=1:nnode 
+		if Connec(iel,i) > 0 loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]]; end
+	end;  
+	Ue=U(loce);
+    sig = feval([deblank(Typel(iel,:)),'_stress'],iel,Ue);
+    % contrainte moyenne de Von Mises sur l'element
+    sVM = sqrt(sig(1)^2+sig(2)^2-sig(1)*sig(2)+3*sig(3)^2);   
+    VM(iel)=sVM ;
+    fprintf('Valeur Moyenne des contraintes dans l''element %3i\n',iel)
+    fprintf('sigxx = %8.3e  sigyy = %8.3e   sigxy = %8.3e  sigVM = %8.3e \n',sig(1),sig(2),sig(3),sVM)
+end
+plot_sig(VM)
+return

Data/NUM_exo2.m

@@ -0,0 +1,128 @@
+close all ; clc;
+%  NUM_exo2
+%  Poutre console avec un mailllage automatique en Q4
+%  H.Oudin 
+
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('Elasticite plane : Poutre console avec un mailllage automatique en Q4');
+disp('=================');
+
+% donnees unites N , mm
+L = 3000 ; h =200 ; e=10;    
+E=210.e+3; nu=0.3;
+poid=200;
+
+% ------ maillage
+nex = input('donner le nombre d''element en x (longueur) ? [5]: ');
+if isempty(nex) nex=5; end  
+ney = input('donner le nombre d''element en y (hauteur) ? [3]: ');
+if isempty(ney) ney=3; end           
+Coord=[];                                                
+for i=0:nex
+    for j=0:ney Coord=[Coord;[i*L/nex  j*h/ney]]; end 
+ end
+Connec=[]; pas=ney+1;   
+for i=1:nex
+    for j=1:ney 
+    i1=(i-1)*pas+j;
+    Connec=[Connec;[i1  i1+pas i1+pas+1 i1+1]];
+    end 
+ end
+[nnod,ndim]=size(Coord);
+[nelt,nnode]=size(Connec);
+nddln=2;  nddlt=nddln*nnod;               
+for i=1:nelt  Typel = char('Q4_ep',Typel); end
+
+% ------ caracteristiques mecaniques (E nu e fx fy)
+Nprop=ones(nex*ney);   
+Prop=[ E nu e 0 0 ];
+
+% ------ definition des CL en deplacement
+CL=[];          
+    for i=0:ney  CL=[CL;[i+1,  1 , 1 ]];  end
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);        
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+
+% ------ definition des charges nodales
+Charg=[(ney+1)*(nex+1) , 0, poid ];     
+F=zeros(nddlt,1);	 
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+ end
+% ----- trace du maillage
+%plotstr  
+% ----- resolution du probleme
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   
+plotdef(U)
+VM = [];
+% -----	post-traitement
+
+sige=[];   % tableau des contraintes sigxx
+for iel=1:nelt          %----- boucle sur les elements
+    loce=[]; 
+    for i=1:nnode 
+		if Connec(iel,i) > 0 loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]]; end
+	end;  
+	Ue=U(loce);
+    sig = feval([deblank(Typel(iel,:)),'_stress'],iel,Ue);
+    % contrainte moyenne de Von Mises sur l'element
+    sVM = sqrt(sig(1)^2+sig(2)^2-sig(1)*sig(2)+3*sig(3)^2);   
+    VM(iel)=sVM ;
+    %fprintf('Valeur Moyenne des contraintes dans l''element %3i\n',iel)
+    %fprintf('sigxx = %8.3e  sigyy = %8.3e   sigxy = %8.3e  sigVM = %8.3e \n',sig(1),sig(2),sig(3),sVM)
+    sige=[sige;[sig(1)]];
+end
+disp(' ');disp('------- Contraintes moyenne de Von Mises sur les elements ----------');
+for iel=1:nelt         
+ fprintf('element %3i sigma VM : %8.3e \n',iel,VM(iel));
+end
+taille = get(0,'ScreenSize'); 
+figure('Name','Contrainte moyenne de Von Mises sur les elements',...
+        'Position',[taille(3)/2.01 taille(4)/2.8 taille(3)/2 taille(4)/2])  
+plot_sig(VM) 
+
+% comparaison avec la solution poutre
+I=e*h^3/12; EI=E*I ;        % solution poutre
+v=poid*(L^3)/(3*EI); sigxx=poid*L*h/(2*I);
+disp(' ');
+fprintf('Solution analytique poutre\n')
+fprintf('la fleche en bout est = %8.3e mm et la contrainte max sigxx = %8.3e  MPa \n',v,sigxx)
+disp(' ');
+T=linspace(0,L,1000);
+for i=1:size(T,2)
+    Sol(i)=poid*(T(i)^2*(3*L-T(i)))/(6*EI);
+end
+taille = get(0,'ScreenSize');
+figure('Name','comparaison des fleches 2D-Q4 (rouge) / poutre (bleu)',...
+       'Position',[taille(3)/2.02 taille(4)/2.6 taille(3)/2 taille(4)/2] )  
+plot(T,Sol);
+hold on
+Uv=U([2:2*(ney+1):2*(nex+1)*(ney+1)]);
+Tv=linspace(0,L,size(Uv,1));
+plot(Tv,Uv,'r');
+figure('Name','Contraintes max sigxx 2D-Q4 (rouge) / poutre (bleu)',...
+       'Position',[taille(3)/2.02 taille(4)/2.6 taille(3)/2 taille(4)/2] )   
+hold on,
+% La contrainte sigxx poutre est calculée au centre des éléments de la rangée du bas
+x=0:L; y=poid*h*(1-1/ney)*(L-x)/(2*I); plot(x,y,'b'),
+for i=1:nex      %contrainte sigxx sur les éléments de la rangée du bas
+    x1=(i-1)*L/nex ; x2 = i*L/nex ;  
+    sign=sige(1+ney*(i-1));
+    x=x1:(x2-x1)/20: x2;
+    y=sign*ones(1,length(x));
+    plot(x,y,'r')
+end
+grid 
+return

Data/colonne.m

@@ -0,0 +1,121 @@
+close all ; clc;
+%  Donnees de la colonne etudiee dans les exercices de cours
+%  Colonne
+%  avec h=100 , E = 1000 , S=1 , mg = 6, maillage : 3 elements de degre 1
+% 
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec ne noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   barre_stress    : calcul de la contrainte dans un element barre
+%   barre_compar     : comparaison avec la solution analytique
+%
+% Initialisation des variables globales
+%   nddln : nb de ddl par noeud
+% 	nnod  : nb de noeuds 
+%  	nddlt : nb de ddl total(=ndln*nnod)   
+%	  nelt  : nb d'elements  
+% 	nnode : nb de noeuds par element (2)
+% 	ndim  : dimension du probleme (1D,2D ou 3D)
+%   ncld  : nb de conditions de champ impose (dirichlet)
+%
+% 	Coord(nnod,ndim): coordonnees des noeuds
+%   Connec(nelt,2)	: connectivites	des elements
+%   Typel(nelt)     : Type des elements (barre_ke)
+% 	Nprop(nelt)		: Ne de caracteristique pour chaque element
+%   Prop(nprop,ncar): Tableau des caracteristiques mecaniques (ES, f)     
+%	  Ncl(1,nddlt)	: vaut 1 si le ddl est impose (deplacements imposes)
+%	  Vcl(1,nddlt)	: valeur du deplacement impose 
+%	  F (nddlt,1)		: vecteur des charges nodales donnees
+
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('Structure etudiee : Colonne traitee en exercice de cours N 8');
+disp('==================');
+% definition du maillage
+h = 100;                    % definition des coordonnees des noeuds
+Coord=[0 ; h ; 3*h ; 6*h ];
+[nnod,ndim]=size(Coord);
+nddln=1;  nddlt=nddln*nnod;     
+ 
+Connec=[ 1 , 2 ; ...        % definition de la matrice de connectivite
+         2 , 3 ; ...
+         3 , 4 ];
+[nelt,nnode]=size(Connec);
+
+% definition du modele EF : type des elements
+Typel = 'barre_ke';               % definition du type des elements
+for i=1:nelt
+    Typel = char('barre_ke',Typel);
+end
+
+% definition des caracteristiques elementaires (ES fx)
+Nprop=[1;1;1];              % pour chaque element Ne de la propriete
+Prop =[1000 -0.01];         % tableau des valeurs de ES fx     
+
+% definition des CL en deplacement
+CL=[ 1 , 1 ]; ...        % Ne du noeud, type   (1 ddl impose, 0 ddl libre)
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);      % Valeurs des deplacements imposes
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1; ncld=ncld+1; end
+   end
+end
+
+% definition des charges nodales
+Charg=[ 4  0 ];          %  Ne du noeud , Fx      
+F=zeros(nddlt,1);	     %  vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+
+% trace du maillage pour validation des donnees 
+plotstr   
+reponse = '';                  
+%%reponse = input('Voulez-vous continuer? O/N [O]: ','s');
+if isempty(reponse) | reponse =='O'
+    %clc; echo on 
+    % echo off
+    U = zeros(nddlt,1);
+    R = zeros(nddlt,1);
+    [U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+
+                                %----- format d'impression des vecteurs
+    form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+    disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+    fprintf(format,U)
+    %plotdef(U)   
+                                            %-----	post-traitement
+    disp(' ');disp('------- Efforts aux appuis  ----------');
+    fprintf(format,R(:,1));
+    [Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+
+    disp(' ');
+    fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+    fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+    fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+
+    disp(' ');disp('------- Contraintes sur les elements ----------');
+    for iel=1:nelt          %----- boucle sur les elements
+        loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+        Ue=U(loce);
+        Ne=feval('barre_stress',iel,Ue);
+        fprintf('Dans l''element %3i l''effort normal est %8.3e\n',iel,Ne)
+    end
+
+   reponse = input('Voulez-vous comparer avec la solution analytique? O/N [O]: ','s');
+   if isempty(reponse) | reponse =='O'
+    feval('barre_compar',U);   %----- comparaison avec la solution analytique
+   end                          
+return
+else
+    disp(' ');disp('---------------- arret du calcul----------------');
+    clear all
+end

Data/portique_MEF2.m

@@ -0,0 +1,83 @@
+close all ; clc;
+%  script du portique de l'exercice MEF2
+%                 
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec numero noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   poutre_stress   : calcul de la contrainte dans un element barre
+%  
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F 
+disp(' ');
+disp('structure etudiee : portique de l''exercice MEF2');
+disp('==================');
+% definition du maillage
+h = 1; L = 1;
+Coord=[ 0 , 0 ; ...         % definition des coordonnees des noeuds X , Y
+        0 , h ; ...
+        L , h ; ...
+        L , 0];
+[nnod,ndim]=size(Coord);
+nddln=3;  nddlt=nddln*nnod;     
+Connec=[ 1 , 2 ; ...        % definition de la matrice de connectivite i , j
+         2 , 3 ; ...
+         3 , 4];
+[nelt,nnode]=size(Connec);
+
+% definition du modele EF : type des elements
+Typel = 'poutre_ke';     
+for i=1:nelt Typel = char('poutre_ke',Typel); end
+
+% definition des caracteristiques mecaniques elementaires (ES EI fx fy)
+Nprop=[1;2;1];          
+Prop=[ 1000 1 0 0 ; 1000 1 0 0 ]; 
+     
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 , 1; ...        % 1 ddl impose ,0 ddl libre
+     4 , 1 , 1 , 1];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);         % Valeurs imposees nulles
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+
+% definition des charges nodales
+Charg=[ 2  84  0  0 ];      %  numero du noeud , Fx,Fy,Mz      
+F=zeros(nddlt,1);	          %  vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+ end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+
+plotstr  % trace du maillage pour validation des donnees 
+
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+%-----	post-traitement
+plotdef(U)
+%----- format d'impression des vecteurs
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)         
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt   
+    loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+    Ue=U(loce);
+    feval('poutre_stress',iel,Ue);
+end                     
+return

Data/portique_cours.m

@@ -0,0 +1,82 @@
+close all ; clc;
+%  jeu de donnees du portique traite en exemple dans le cours
+%  avec h=1  , ES = 1000, EI = 1 , F = 30.
+% 
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec numero noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   poutre_stress    : calcul de la contrainte dans un element poutre
+%
+% Initialisation des variables globales pour un portique 2D par H.Oudin
+%   
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F 
+disp(' ');
+disp('structure etudiee : portique traite en exemple dans le cours');
+disp('==================');
+% definition du maillage
+h = 1;
+Coord=[ 0 , 0 ; ...         % definition des coordonnees des noeuds X , Y
+        0 , h ; ...
+        h , h ];
+[nnod,ndim]=size(Coord);
+nddln=3;  nddlt=nddln*nnod;    
+
+Connec=[ 1 , 2 ; ...        % definition de la matrice de connectivite i , j
+         2 , 3 ];
+[nelt,nnode]=size(Connec);
+
+% definition du modele EF : type des elements
+Typel = 'poutre_ke'; 
+for i=1:nelt
+  Typel = char('poutre_ke',Typel);
+end    
+% definition des caracteristiques mecaniques elementaires (ES fx fy)
+Nprop=[1;1];              % pour chaque element numero de la propriete
+Prop=[ 1000 1 0 0 ];      % tableau des differentes valeurs de ES EI fx fy    
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 , 1; ...   % numero du noeud, type sur u,v,teta (1 ddl impose ,0 ddl libre)
+     3 , 0 , 1 , 1];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);       % Valeurs imposees nulles
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+% definition des charges nodales
+Charg=[ 2  30  0  0  ];    %  numero du noeud , Fx,Fy,Mz  
+F=zeros(nddlt,1);	         %  vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+ end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+
+plotstr  % trace du maillage pour validation des donnees 
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)
+plotdef(U)                     %-----	post-traitement
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+   loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+   Ue=U(loce);
+   feval('poutre_stress',iel,Ue);
+end                       
+return

Data/portique_exo17.m

@@ -0,0 +1,80 @@
+close all ; clc;
+%  Script pour la portique de l'exercice de cours exo17
+%  maillage  avec nelt elements
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec numero noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   poutre_stress   : calcul de la contrainte dans un element barre
+%   
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F 
+disp(' ');
+disp('structure etudiee : portique de l''exercice de cours exo17');
+disp('==================');
+disp('les 2 poutres sont maillees en ne elements');
+ne = input('donner le nombre d''elements ne ? [1]: ');
+if isempty(ne) ne=1; end 
+% definition du maillage
+h = 1; nelt=2*ne;
+Coord=[]; 
+for j=0:ne Coord=[Coord;[0  j*h/ne]]; end 
+for j=1:ne Coord=[Coord;[j*h/ne h ]]; end 
+[nnod,ndim]=size(Coord);
+nddln=3;  nddlt=nddln*nnod;  
+Connec=[]; nnode = 2;
+for j=1:nelt Connec=[Connec;[j  j+1]]; end    
+% definition du modele EF : type des elements
+Typel = 'poutre_ke';         
+for i=1:nelt Typel = char('poutre_ke',Typel); end
+% definition des caracteristiques mecaniques elementaires (ES EI fx fy)
+Nprop=[]; 
+for j=1:ne Nprop=[Nprop; 1]; end 
+for j=1:ne Nprop=[Nprop; 2]; end 
+Prop=[ 2 1 -560 0          % tableau des differentes valeurs de ES EI fx fy
+       2 1 0   0 ];     
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 , 1; ...    % numero du noeud, type (1 ddl impose ,0 ddl libre)
+    nnod , 1 , 1 , 1];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);         % deplacements imposes nuls
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+% definition des charges nodales
+Charg=[ 1  0  0  0               %  numero du noeud , Fx,Fy,Mz
+      ];
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+ end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+plotstr  % trace du maillage
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+%-----	post-traitement
+plotdef(U)
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)                                       
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+    loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+    Ue=U(loce);
+    feval('poutre_stress',iel,Ue);
+end                     
+return

Data/portique_exo17b.m

@@ -0,0 +1,79 @@
+close all ; clc;
+%  Script pour la portique de l'exercice de cours exo17
+%  maillage  avec nelt elements
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec numero noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   poutre_stress   : calcul de la contrainte dans un element barre
+%   
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F 
+disp(' ');
+disp('structure etudiee : portique de l''exercice de cours exo17');
+disp('==================');
+disp('la poutre verticale est maillee en ne elements');
+ne = input('donner le nombre d''elements ne ? [1]: ');
+if isempty(ne) ne=1; end 
+% definition du maillage
+h = 1; nelt=ne+1;
+Coord=[]; 
+for j=0:ne Coord=[Coord;[0  j*h/ne]]; end 
+Coord=[Coord;[h  h]];
+[nnod,ndim]=size(Coord);
+nddln=3;  nddlt=nddln*nnod;  
+Connec=[]; nnode = 2;
+for j=1:nelt Connec=[Connec;[j  j+1]]; end    
+% definition du modele EF : type des elements
+Typel = 'poutre_ke';         
+for i=1:nelt Typel = char('poutre_ke',Typel); end
+% definition des caracteristiques mecaniques elementaires (ES EI fx fy)
+Nprop = ones(nelt);
+Nprop(ne+1)=2;  
+Prop=[ 2 1 -560 0          % tableau des differentes valeurs de ES EI fx fy
+       2 1 0   0 ];     
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 , 1; ...    % numero du noeud, type (1 ddl impose ,0 ddl libre)
+     nnod , 1 , 1 , 1];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);         % deplacements imposes nuls
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+% definition des charges nodales
+Charg=[ 1  0  0  0               %  numero du noeud , Fx,Fy,Mz
+      ];
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+ end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+plotstr  % trace du maillage
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+%-----	post-traitement
+plotdef(U)
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)                                       
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+    loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+    Ue=U(loce);
+    feval('poutre_stress',iel,Ue);
+end                     
+return

Data/poutre_exo15.m

@@ -0,0 +1,76 @@
+close all ; clc;
+%  Script pour la poutre de l'exercice de cours exo15
+%         poutre sur 2 appuis, modelisee par nelt elements
+% 
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec numero noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   poutre_stress   : calcul de la contrainte dans un element poutre
+%   poutre_compar   : comparaison avec la solution analytique
+%   
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F 
+disp(' ');
+disp('structure etudiee : poutre de l''exercice de cours exo15');
+disp('==================');
+% definition du maillage
+nelt = input('donner le nombre d''elements ne ? [2]: ');
+if isempty(nelt) nelt=2; end 
+L = 1;
+Coord=[]; 
+for j=0:nelt Coord=[Coord; j*L/nelt]; end 
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;
+Connec=[]; nnode = 2;
+for j=1:nelt Connec=[Connec;[j  j+1]]; end       
+% definition du modele EF : type des elements
+Typel = 'poutre_ke';       
+for i=1:nelt Typel = char('poutre_ke',Typel); end
+% definition des caracteristiques mecaniques elementaires (EI f)  (en 1D)
+Nprop = ones(nelt);   % pour chaque element numero de la propriete
+Prop=[ 1 -384 ];         % tableau des differentes valeurs de EI fy    
+% definition des CL en deplacement
+CL=[ 1 , 1 , 0 ; ...     % numero du noeud, (1 ddl impose ,0 ddl libre)
+    nelt+1 , 1 , 0 ];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);         % Valeurs imposees nulles
+for i=1:size(CL,1)
+    for j=1:nddln 
+    if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+    end
+end
+% definition des charges nodales
+F=zeros(nddlt,1);	   
+[Fx,Fy,Fz] = feval('resultante',F); 
+
+plotstr  % trace du maillage pour validation des donnees
+% ----- resolution du probleme
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   
+plotdef(U)
+%-----	post-traitement
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)                                     
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+  loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+  Ue=U(loce);
+  feval('poutre_stress',iel,Ue);
+  end
+reponse = input('Voulez-vous comparer avec la solution analytique? O/N [O]: ','s');
+if isempty(reponse) | reponse =='O'
+ feval('poutre_compar',U); %----- comparaison avec la solution analytique
+ end                                               
+clear all
+return

Data/treillis_MEF6.m

@@ -0,0 +1,81 @@
+close all ; clc;
+%  Treillis de l'exercice MEF6
+%  H.Oudin 
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('Structure etudiee  : Treillis de l''exercice MEF6 ');
+disp('=================');
+% definition du maillage
+h = 100*sqrt(2)/2;   % definition des coordonnees des noeuds X , Y
+Coord=[ 0   0         
+   h  -h 
+   2*h  0 
+   3*h  -h
+   4*h  0 ];
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;     
+% definition de la matrice de connectivite i , j
+Connec=[ 1 2 ; 1 3 ; 2 3 ; 2 4 ; 3 4 ; 3 5 ; 4 5 ];
+[nelt,nnode]=size(Connec); 
+
+Typel = 'barre_ke';               % definition du type des elements
+for i=1:nelt
+    Typel = char('barre_ke',Typel);
+end
+
+Nprop=ones(nelt);
+Prop=[ 100*sqrt(2) 0 0 ];       % tableau des differentes valeurs de ES fx fy      
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 ; ...      
+     5 , 0 , 1 ];
+Ncl=zeros(1,nddlt); ncld=0;
+Vcl=zeros(1,nddlt);       % Valeurs des deplacements imposes
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 
+           Ncl(1,(CL(i,1)-1)*nddln+j)=1;
+           ncld=ncld+1; 
+       end
+   end
+end
+
+% definition des charges nodales    numero du noeud , Fx , Fy
+Charg=[ 2    0  -1.                
+        4    0  -1.];
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+     F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+  end
+end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+% trace du maillage pour validation des donnees 
+plotstr                     
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+
+%----- format d'impression des vecteurs
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)
+plotdef(U)
+%-----	post-traitement
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+   loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+   Ue=U(loce);
+   Ne=feval('barre_stress',iel,Ue);
+   fprintf('Dans l''element %3i l''effort normal est %8.3e\n',iel,Ne)
+   end  
+return

Data/treillis_MEF6_sym.m

@@ -0,0 +1,78 @@
+close all ; clc;
+%  Treillis de l'exercice MEF6 en tenant compte de la symetrie
+%  H.Oudin 
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('Structure etudiee  : Exercice MEF6 avec symetrie');
+disp('=================');
+% definition du maillage
+h = 100*sqrt(2)/2;  
+Coord=[ 0   0; h  -h ; 2*h  0 ; 2*h  -h ];
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;     
+Connec=[ 1 2 ; 1 3 ; 2 3 ; 2 4 ];
+[nelt,nnode]=size(Connec); 
+Typel = 'barre_ke';              
+for i=1:nelt
+    Typel = char('barre_ke',Typel);
+end
+Nprop=ones(nelt);
+Prop=[100*sqrt(2) 0 0 ];   
+CL=[ 1 , 0 , 1 ; ...     
+     3 , 1 , 0 ; ...
+     4 , 1 , 1 ; ...
+     ]; 
+Ncl=zeros(1,nddlt); ncld=0;
+Vcl=zeros(1,nddlt);     
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 
+           Ncl(1,(CL(i,1)-1)*nddln+j)=1;
+           ncld=ncld+1; 
+       end
+   end
+end
+Charg=[ 2    0  -1.];
+F=zeros(nddlt,1);	    
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+end
+[Fx,Fy,Fz] = feval('resultante',F);     
+
+% trace du maillage 
+plotstr                     
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+
+%----- format d'impression des vecteurs
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)
+plotdef(U)
+%-----	post-traitement
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+   loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+   Ue=U(loce);
+   Ne=feval('barre_stress',iel,Ue);
+   fprintf('Dans l''element %3i l''effort normal est %8.3e\n',iel,Ne)
+   end  
+disp(' ');   
+disp('Expliquer pourquoi les resultats sont justes,');
+disp('alors que le noeud 4 est bloque !!');
+return
+
+

Data/treillis_cours.m

@@ -0,0 +1,122 @@
+close all ; clc;
+%  Script de calcul statique du treillis traite en exemple dans le cours
+%  avec a=100 , ES = 100*sqrt(2) , F = 40.
+% 
+% Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec ne noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   barre_stress    : calcul de la contrainte dans un element barre
+%
+% Initialisation des variables globales pour un treillis
+%   nddln : nb de ddl par noeud
+% 	nnod  : nb de noeuds 
+%  	nddlt : nb de ddl total(=ndln*nnod)   
+%	  nelt  : nb d'elements  
+% 	nnode : nb de noeuds par element (2)
+% 	ndim  : dimension du probleme (1D,2D ou 3D)
+%   ncld  : nb de conditions de champ impose (dirichlet)
+%
+% 	Coord(nnod,ndim): coordonnees des noeuds
+%   Connec(nelt,2)	: connectivites	des elements
+%   Typel(nelt)     : Type des elements (barre_ke)
+% 	Nprop(nelt)		  : Ne de caracteristique pour chaque element
+%   Prop(nprop,ncar): Tableau des caracteristiques mecaniques (ES, f)     
+%	  Ncl(1,nddlt)	  : vaut 1 si le ddl est impose (deplacements imposes)
+%	  Vcl(1,nddlt)	  : valeur du deplacement impose 
+%	  F (nddlt,1)		  : vecteur des charges nodales donnees
+%
+% L'objectif de scripts de donnees est d'initialiser ses variables globales,
+% avant de lancer les calculs et d'exploiter les resultats (post-traitement).
+% 
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('structure etudiee : treillis traite en exemple dans le cours');
+disp('==================');
+% definition du maillage
+h = 100*sqrt(2)/2;
+Coord=[ 0 , 0 ; ...         % definition des coordonnees des noeuds X , Y
+        2*h , 0 ; ...
+        h , h ];
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;     
+ 
+Connec=[ 1 , 2 ; ...        % definition de la matrice de connectivite i , j
+         1 , 3 ; ...
+         2 , 3 ];
+[nelt,nnode]=size(Connec);
+
+%definition du modele EF : type des elements
+Typel = 'barre_ke';               % definition du type des elements
+for i=1:nelt
+    Typel = char('barre_ke',Typel);
+end
+% definition des caracteristiques mecaniques elementaires (ES fx fy)
+Nprop=[1;1;1];              % pour chaque element Ne de la propriete
+Prop=[ 100*sqrt(2) 0 0      % tableau des differentes valeurs de ES fx fy
+      ];
+      
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 ; ...  % Ne du noeud, type sur u et v (1 ddl impose ,0 ddl libre)
+     2 , 0 , 1 ];
+Ncl=zeros(1,nddlt); ncld=0;
+Vcl=zeros(1,nddlt);  % Valeurs des deplacements imposes
+%Vcl(2)=1;           % à utiliser pour imposer une valeur non nulle sur un ddl i
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 
+           Ncl(1,(CL(i,1)-1)*nddln+j)=1;
+           ncld=ncld+1; 
+       end
+   end
+end
+
+% definition des charges nodales
+Charg=[ 3  40.  0                 %  Ne du noeud , Fx , Fy
+      ];
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+
+disp('Les variables globales sont initialisees');
+disp('Fin de lecture des donnees');
+% trace du maillage pour validation des donnees 
+plotstr                     
+reponse = input('Voulez-vous continuer? O/N [O]: ','s');
+if isempty(reponse) | reponse =='O'
+    U = zeros(nddlt,1);
+    R = zeros(nddlt,1);
+    [U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+%----- format d'impression des vecteurs
+    form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+    disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+    fprintf(format,U)
+    plotdef(U)
+%-----	post-traitement
+    disp(' ');disp('------- Efforts aux appuis  ----------');
+    fprintf(format,R(:,1));
+    [Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+    disp(' ');
+    fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+    fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+    fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+    disp(' ');disp('------- Contraintes sur les elements ----------');
+    for iel=1:nelt          %----- boucle sur les elements
+        loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+        Ue=U(loce);
+        Ne = feval('barre_stress',iel,Ue);
+        fprintf('Dans l''element %3i l''effort normal est %8.3e\n',iel,Ne)
+    end 
+    clear all
+    return
+else
+disp(' ');disp('---------------- arret du calcul----------------');
+clear all
+end

Data/treillis_cours_V2.m

@@ -0,0 +1,96 @@
+close all ; clc;
+%  jeu de donnees du treillis traite en exemple dans le cours
+%  avec a=100 , ES = 100*sqrt(2) , F = 40.
+%  avec la methode du terme unite sur la diagonale
+%
+% Fonctions utilisees
+%   statique        : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec ne noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   barre_stress    : calcul de la contrainte dans un element barre
+
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+disp(' ');
+disp('structure etudiee : treillis traite en exemple dans le cours');
+disp('==================');
+% definition du maillage
+h = 100*sqrt(2)/2;
+Coord=[ 0 , 0 ; ...         % definition des coordonnees des noeuds X , Y
+        2*h , 0 ; ...
+        h , h ];
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;     
+ 
+Connec=[ 1 , 2 ; ...        % definition de la matrice de connectivite i , j
+         1 , 3 ; ...
+         2 , 3 ];
+[nelt,nnode]=size(Connec);
+
+%definition du modele EF : type des elements
+Typel = 'barre_ke';               % definition du type des elements
+for i=1:nelt
+    Typel = char('barre_ke',Typel);
+end
+% definition des caracteristiques mecaniques elementaires (ES fx fy)
+Nprop=[1;1;1];              % pour chaque element numero de la propriete
+Prop=[ 100*sqrt(2) 0 0      % tableau des differentes valeurs de ES fx fy
+      ];
+      
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 ; ... % numero du noeud, type sur u et v (1 ddl impose ,0 ddl libre)
+     2 , 0 , 1 ];
+Ncl=zeros(1,nddlt); ncld=0;
+Vcl=zeros(1,nddlt);       % Valeurs des deplacements imposes
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 
+           Ncl(1,(CL(i,1)-1)*nddln+j)=1;
+           ncld=ncld+1; 
+       end
+   end
+end
+
+% definition des charges nodales
+Charg=[ 3  40.  0                 %  numero du noeud , Fx , Fy
+      ];
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+    F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+    end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+
+% trace du maillage pour validation des donnees 
+close all
+plotstr                      
+U = zeros(nddlt,1);
+R = zeros(nddlt,1);
+[U(:,1),R(:,1)] = statique;   % ----- resolution du probleme
+%                               methode du terme unite sur la diagonale                              
+%----- format d'impression des vecteurs
+form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+fprintf(format,U)
+plotdef(U)                                             
+%-----	post-traitement
+disp(' ');disp('------- Efforts aux appuis  ----------');
+fprintf(format,R(:,1));
+[Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+
+disp(' ');
+fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+
+disp(' ');disp('------- Contraintes sur les elements ----------');
+for iel=1:nelt          %----- boucle sur les elements
+  loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+  Ue=U(loce);
+  Ne=feval('barre_stress',iel,Ue);
+  fprintf('Dans l''element %3i l''effort normal est %8.3e\n',iel,Ne)
+  end
+clear all

Data/treillis_exo10.m

@@ -0,0 +1,125 @@
+close all ; clc;
+%  treillis : donnees de l'exercice de cours 10 
+%  avec h=1+sqrt(2)  , ES = 100*sqrt(2) , F = 10.
+% 
+%  Fonctions utilisees
+%   statiqueUR      : calcul de la reponse statique [U,R]
+%   plotstr         : trace du maillage avec ne noeuds et elements
+%   plotdef         : trace de la deformee
+%   resultante      : calcul de la resultante en x, y, z d'un vecteur nodal
+%   barre_stress    : calcul de la contrainte dans un element barre
+%
+% Initialisation des variables globales
+%   nddln : nb de ddl par noeud
+% 	nnod  : nb de noeuds 
+%  	nddlt : nb de ddl total(=ndln*nnod)   
+%	  nelt  : nb d'elements  
+% 	nnode : nb de noeuds par element (2)
+% 	ndim  : dimension du probleme (1D,2D ou 3D)
+%   ncld  : nb de conditions de champ impose (dirichlet)
+%
+% 	Coord(nnod,ndim): coordonnees des noeuds
+%   Connec(nelt,2)	: connectivites	des elements
+%   Typel(nelt)     : Type des elements (barre_ke)
+% 	Nprop(nelt)		: Ne de caracteristique pour chaque element
+%   Prop(nprop,ncar): Tableau des caracteristiques mecaniques (ES, f)     
+%	  Ncl(1,nddlt)	: vaut 1 si le ddl est impose (deplacements imposes)
+%	  Vcl(1,nddlt)	: valeur du deplacement impose 
+%	  F (nddlt,1)		: vecteur des charges nodales donnees
+%  H.Oudin 
+  
+global nddln nnod nddlt nelt nnode ndim ncld
+global Coord Connec Typel Nprop Prop Ncl Vcl F 
+
+% definition du maillage
+h = 1+sqrt(2);
+Coord=[ h , -h ; ...         % definition des coordonnees des noeuds X , Y
+        0 , 0 ; ...
+        h , 0 ; ...
+      2*h , 0 ];
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;     
+ 
+Connec=[ 1 , 2 ; ...        % definition de la matrice de connectivite i , j
+         1 , 3 ; ...
+         1 , 4 ];
+[nelt,nnode]=size(Connec);
+
+% definition des caracteristiques mecaniques elementaires (ES fx fy)
+Nprop=[1;1;1];              % pour chaque element Ne de la propriete
+Prop=[ 100*sqrt(2) 0 0      % tableau des differentes valeurs de ES fx fy
+      ];
+
+  % definition du modele EF : type des elements
+Typel = 'barre_ke';               % definition du type des elements
+for i=1:nelt
+    Typel = char('barre_ke',Typel);
+end
+      
+% definition des CL en deplacement
+CL=[ 2 , 1 , 1 ; ...        % Ne du noeud, type sur u et v (1 ddl impose ,0 ddl libre)
+     3 , 1 , 1 ; ...
+     4 , 1 , 1 ];
+Ncl=zeros(1,nddlt); ncld=0;
+Vcl=zeros(1,nddlt);         % Valeurs des deplacements imposes
+%Vcl(2)=1;                  % e utiliser pour imposer une valeur non nulle sur un ddl i
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1;ncld=ncld+1; end
+   end
+end
+
+% definition des charges nodales
+Charg=[ 1  0  -10                 %  Ne du noeud , Fx , Fy
+      ];
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+for iclf=1:size(Charg,1)           
+	noeud=Charg(iclf,1);  
+	for i=1:nddln
+       F((noeud-1)*nddln+i)=F((noeud-1)*nddln+i) + Charg(iclf,i+1);
+    end
+end
+[Fx,Fy,Fz] = feval('resultante',F);      %----- resultante des charges nodales
+clc;disp(' ');
+disp('structure etudiee : treillis de l''exercice de cours N 10');
+disp('==================');
+disp('Les variables globales sont initialisees');
+disp('Fin de lecture des donnees');
+
+% trace du maillage pour validation des donnees 
+plotstr                     
+reponse = input('Voulez-vous continuer? O/N [O]: ','s');
+if isempty(reponse) | reponse =='O'
+    disp(' ');
+    U = zeros(nddlt,1);
+    R = zeros(nddlt,1);
+    [U(:,1),R(:,1)] = statiqueUR;   % ----- resolution du probleme
+
+                                %----- format d'impression des vecteurs
+    form =' %8.3e   %8.3e   %8.3e  '; format = [form(1:8*nddln),' \n']; 
+    disp(' ');disp('------- deplacements nodaux sur (x,y,z) ----------');
+    fprintf(format,U)
+    plotdef(U)
+                                            %-----	post-traitement
+    disp(' ');disp('------- Efforts aux appuis  ----------');
+    fprintf(format,R(:,1));
+    [Rx,Ry,Rz] = feval('resultante',R);     %----- resultantes et reactions
+
+    disp(' ');
+    fprintf('La resultante des charges nodales    en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Fx,Fy,Fz);                    
+    fprintf('La resultante des charges reparties  en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',-Rx-Fx,-Ry-Fy,-Rz-Fz);
+    fprintf('La resultante des efforts aux appuis en (x,y,z) est : %8.3e   %8.3e   %8.3e \n',Rx,Ry,Rz);
+
+    disp(' ');disp('------- Contraintes sur les elements ----------');
+    for iel=1:nelt          %----- boucle sur les elements
+        loce=[]; for i=1:nnode loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];end
+        Ue=U(loce);
+        Ne=feval('barre_stress',iel,Ue);
+        fprintf('Dans l''element %3i l''effort normal est %8.3e\n',iel,Ne)
+    end                       
+clear all
+return
+else
+    disp(' ');disp('---------------- arret du calcul----------------');
+    clear all
+end

Data/vibra_barre_coursMEF.m

@@ -0,0 +1,81 @@
+close all ; clear; clc;
+%  Script de calcul des modes de vibration de la barre encastree - libre
+%  exemple du cours modelisee par nelt elements
+% 
+%  fonction utilisee
+%    plotstr  : trace du maillage avec numero noeuds et elements
+%    vibrations : resolution du probleme
+%
+%  H.Oudin 
+disp(' ');
+disp('modes de vibration d''une barre encastree - libre');
+disp('==================');
+global nddln nnod nddlt nelt nnode ndim
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+
+% definition du maillage
+nelt = input('donner le nombre d''elements nelt ? [2]: ');
+if isempty(nelt) nelt=2; end 
+L = 1;
+Coord=[]; 
+for j=0:nelt Coord=[Coord; j*L/nelt]; end 
+[nnod,ndim]=size(Coord);
+nddln=1;  nddlt=nddln*nnod;
+Connec=[]; nnode = 2;
+for j=1:nelt Connec=[Connec;[j  j+1]]; end 
+
+% definition du modele EF : type des elements
+Typel = 'barre_keme';         
+for i=1:nelt Typel = char('barre_keme',Typel); end
+
+% definition des caracteristiques mecaniques  (en 1D)
+Nprop = ones(nelt);   % pour chaque element numero de la propriete
+Prop=[ 1 1 ]; % tableau de ES RhoS   
+
+% definition des CL en deplacement
+CL=[ 1 , 1 ];
+Ncl=zeros(1,nddlt);
+Vcl=zeros(1,nddlt);         % deplacements imposes nuls
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1; end
+   end
+end
+F=zeros(nddlt,1);
+%plotstr  % trace du maillage                    
+nmode = input('combien de frequences faut-il calculer ? [2]: ');
+if isempty(nmode) nmode=2; end   
+f=[];  U=[]; 
+[f,U] = feval('vibrations',nmode);
+disp('frequences obtenues par le modele EF');
+f
+if f(1)>f(2) 
+  ft=fliplr(f');
+  U=fliplr(U);
+  f=ft';
+end
+
+ES=Prop(1,1);RS=Prop(1,2);
+norme = RS*L/2; 
+taille = get(0,'ScreenSize'); 
+figure('Name','Comparaison analytique / numerique des modes de vibration de la poutre',...
+        'Position',[taille(3)/2.01 taille(4)/2.4 taille(3)/2 taille(4)/2])
+
+for j= 1: nmode 
+   x=[0:0.01*L:L];  y = sin((2*j-1)*pi*x/(2*L));%----- sol analytique appuyeee - appuyee
+   fanal = (2*j-1)*sqrt(ES/(RS*(L^2)))/4;err=100*(f(j)-fanal)/fanal;
+   subplot(nmode,1,j),title([int2str(j),' mode de vibration a la frequence ',...
+   num2str(f(j),'%7.2f'),'  / sol analytique ', num2str(fanal,'%7.2f'), ...
+   '  soit une erreur de ', num2str(err,'%4.2f'), ' %']),hold on
+   plot(x,y,'g'); plot(x,-y,'g'); axis([0  L  -1.5  1.5])
+   s = sqrt(norme);%---- facteur d'echelle sur la deformee EF a partir de la norme analytique
+   for iel = 1:nelt
+    Pos  = Coord(Connec(iel,:),:); Le=Pos(2,1)-Pos(1,1);
+    je=(Connec(iel,1)-1)*nddln+1;
+    x = [0:0.1:1];
+    y = s*( (1-x)*U(je,j) + x.*U(je+1,j));
+    x = Pos(1,1)+ x*Le;
+    plot(x,y,'b'); plot(x,-y,'b')
+   end       
+end
+return

Data/vibra_portique_exo17.m

@@ -0,0 +1,52 @@
+close all ; clc;
+%  Script de calcul des modes de vibration du portique de l'exo 17
+%  portique 2D  avec 2*ne elements
+% 
+%  fonction utilisee
+%    plotstr  : trace du maillage avec numero noeuds et elements
+%    vibrations : resolution du probleme
+%    plotmodes : trace des modes de vibrations d'un portique 2D
+%
+%  H.Oudin 
+
+global nddln nnod nddlt nelt nnode ndim
+global Coord Connec Nprop Prop Ncl Vcl F
+disp(' ');
+disp('modes de vibration du portique de l''exercice de cours 17');
+disp('==================');
+disp('les poutres sont maillees en ne elements');
+ne = input('donner le nombre d''elements ne ? [2]: ');
+if isempty(ne) ne=2; end 
+% definition du maillage
+h = 1; nelt=2*ne;
+Coord=[]; 
+for j=0:ne Coord=[Coord;[0  j*h/ne]]; end 
+for j=1:ne Coord=[Coord;[j*h/ne h ]]; end 
+[nnod,ndim]=size(Coord);
+nddln=3;  nddlt=nddln*nnod;  
+Connec=[]; nnode = 2;
+for j=1:nelt Connec=[Connec;[j  j+1]]; end    
+% tableau des proprietees ES EI rhoS
+Nprop = ones(nelt);
+Prop=[ 21000000 210 0.78 ];      
+% definition des CL en deplacement
+CL=[ 1 , 1 , 1 , 1; ...    % numero du noeud, type (1 ddl impose ,0 ddl libre)
+     nnod , 1 , 1 , 1];
+Ncl=zeros(1,nddlt);ncld=0;
+Vcl=zeros(1,nddlt);         % deplacements imposes nuls
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1; end
+   end
+end
+F=zeros(nddlt,1);	    %----- vecteur sollicitation
+%plotstr  % trace du maillage                    
+nmode = input('combien de frequences faut-il calculer ? [3]: ');
+if isempty(nmode) nmode=3; end   
+[f,U] = feval('vibrations',nmode);
+for j= 1: nmode 
+   figure('Name','modes de vibration du portique',...
+        'Position',[taille(3)/2.01 taille(4)/2.4 taille(3)/2 taille(4)/2])
+   plotmodes(U(:,nmode+1-j),j,f(nmode+1-j))
+end
+return

Data/vibra_poutre_exo15.m

@@ -0,0 +1,84 @@
+close all ; clc;
+%  Script de calcul des modes de vibration de la poutre de l'exo 15
+%  exemple du cours modelisee par nelt elements
+% 
+%  fonction utilisee
+%    plotstr  : trace du maillage avec numero noeuds et elements
+%    vibrations : resolution du probleme
+%
+%  H.Oudin 
+disp(' ');
+disp('modes de vibration d''une poutre appuyee - appuyee');
+disp('==================');
+global nddln nnod nddlt nelt nnode ndim
+global Coord Connec Typel Nprop Prop Ncl Vcl F
+
+% definition du maillage
+nelt = input('donner le nombre d''elements nelt ? [2]: ');
+if isempty(nelt) nelt=2; end 
+L = 1;
+Coord=[]; 
+for j=0:nelt Coord=[Coord; j*L/nelt]; end 
+[nnod,ndim]=size(Coord);
+nddln=2;  nddlt=nddln*nnod;
+Connec=[]; nnode = 2;
+for j=1:nelt Connec=[Connec;[j  j+1]]; end 
+
+% definition du modele EF : type des elements
+Typel = 'poutre_keme';         
+for i=1:nelt Typel = char('poutre_keme',Typel); end
+
+% definition des caracteristiques mecaniques elementaires (EI f)  (en 1D)
+Nprop = ones(nelt);   % pour chaque element numero de la propriete
+Prop=[ 21000000 210 0.78 ]; % tableau de ES EI RhoS   
+
+% definition des CL en deplacement
+CL=[ 1 , 1 , 0 ; ...     % numero du noeud, (1 ddl impose ,0 ddl libre)
+    nelt+1 , 1 , 0 ];
+Ncl=zeros(1,nddlt);
+Vcl=zeros(1,nddlt);         % deplacements imposes nuls
+for i=1:size(CL,1)
+   for j=1:nddln 
+       if CL(i,1+j)==1 Ncl(1,(CL(i,1)-1)*nddln+j)=1; end
+   end
+end
+F=zeros(nddlt,1);
+%plotstr  % trace du maillage                    
+nmode = input('combien de frequences faut-il calculer ? [3]: ');
+if isempty(nmode) nmode=3; end
+f=[];  U=[]; 
+[f,U] = feval('vibrations',nmode);
+disp('frequences obtenues par le modele EF');
+f
+if f(1)>f(2) 
+  ft=fliplr(f');
+  U=fliplr(U);
+  f=ft';
+  end
+
+EI=Prop(1,2);RS=Prop(1,3);
+norme = RS*L/2;
+disp('en bleu les modes EF / en vert la solution analytique');
+taille = get(0,'ScreenSize'); 
+figure('Name','Comparaison analytique / numerique des modes de vibration de la poutre',...
+        'Position',[taille(3)/2.01 taille(4)/2.4 taille(3)/2 taille(4)/2])
+
+for j= 1: nmode 
+   x=[0:0.01*L:L];  y = sin(j*pi*x/L);%----- sol analytique appuyeee - appuyee
+   fanal = (j^2)*sqrt(EI/(RS*(L^4)))*pi/2;
+   subplot(nmode,1,j),title([int2str(j),' mode de vibration a la frequence ',...
+   num2str(f(j),'%7.2f'),'  / sol analytique ', ...
+   num2str(fanal,'%7.2f')]),hold on
+   plot(x,y,'g'); plot(x,-y,'g'); axis([0  L  -1.2  1.2])
+   s = sqrt(norme);%---- facteur d'echelle sur la deformee EF a partir de la norme analytique
+   for iel = 1:nelt
+    Pos  = Coord(Connec(iel,:),:); Le=Pos(2,1)-Pos(1,1);
+    je=(Connec(iel,1)-1)*nddln+1;
+    x = [0:0.1:1]; x2=x.*x ; x3=x.*x2;
+    y = s*( (1-3*x2+2*x3)*U(je,j) + Le*(x-2*x2+x3)*U(je+1,j) + ...
+    (3*x2-2*x3)*U(je+2,j) + Le*(-x2+x3)*U(je+3,j) );
+    x = Pos(1,1)+ x*Le;
+    plot(x,y,'b'); plot(x,-y,'b')
+   end       
+end
+return

Dessin/plot2D.m

@@ -0,0 +1,50 @@
+function plot2D
+% trace du maillage d'une structure plane 
+% avec les numero des noeuds et des elements
+% les conditions aux limites (deplacements imposes) 
+% et les Charges nodales
+%
+% appel plot2D
+%
+%  H.Oudin 
+global nddln nnod nddlt nelt nnode ndim
+global Coord Connec Ncl F
+
+hold on
+title('maillage de la structure')
+axis equal
+
+XY = Coord;
+Cu = Ncl(1:nddln:nddlt); Cv = Ncl(2:nddln:nddlt);
+Fx = F(1:nddln:nddlt);   Fy = F(2:nddln:nddlt);
+if ndim == 1 
+    XY(:,2) = zeros(nnod,1); 
+    if (nddln == 1) Cv = zeros(nnod,1); Fy = zeros(nnod,1); end
+end
+  
+for inod=1:nnod     % ----- Noeuds , conditions aux limites et charges nodales
+    text(XY(inod,1),XY(inod,2),[' ',int2str(inod)],'color','m','FontSize',14)
+    if Cu(inod)==1 plot(XY(inod,1),XY(inod,2),'>','color','k','MarkerSize',10);end
+    if Cv(inod)==1 plot(XY(inod,1),XY(inod,2),'^','color','k','MarkerSize',10);end
+    if Fx(inod)~= 0  
+        text(XY(inod,1),XY(inod,2),['  \rightarrow',num2str(Fx(inod))],...
+            'HorizontalAlignment','left','color','g','FontSize',14);
+    end
+    if Fy(inod)~= 0  
+        text(XY(inod,1),XY(inod,2),['\uparrow',num2str(Fy(inod))],...
+            'VerticalAlignment','bottom','color','g','FontSize',14);
+    end   
+end
+plot(XY(:,1),XY(:,2),'.','color','m','MarkerSize',12);
+
+for iel = 1:nelt        %----- visualisation du maillage
+    loce=[];            %----- table de localisation pour l'element
+    for i=1:nnode 
+		if Connec(iel,i) > 0 loce=[loce,Connec(iel,i)]; end
+	end;  
+    Pos  = XY(loce,:);
+    X = [[Pos(:,1)]; Pos(1,1)];Y = [[Pos(:,2)]; Pos(1,2)];
+    line(X,Y,'color','b')
+    text(mean(Pos(:,1)),mean(Pos(:,2)),[' ',int2str(iel)],'color','b','FontSize',14)
+end 
+return

Dessin/plotdef.m

@@ -0,0 +1,56 @@
+function  plotdef(Sol)
+% dessin de la deformee d'une structure / a la position initiale 
+%
+% appel plotdef(Sol)
+%       en entree Sol : vecteur des deplacements nodaux dimension (nddlt)
+%
+%  H.Oudin
+%==========================================================================
+global nddln  nnod nddlt nelt nnode ndim
+global Coord Connec
+taille = get(0,'ScreenSize'); 
+figure('Name','deformee de la structure',...
+        'Position',[taille(3)/2.02 taille(4)/2.6 taille(3)/2 taille(4)/2])
+axis equal
+    grid on
+    xlabel('X');
+    ylabel('Y');
+    zlabel('Z');    
+
+XY = Coord;
+if ndim == 1 
+    XY(:,2) = zeros(nnod,1);XY(:,3) = zeros(nnod,1); 
+    if nddln == 2 
+        u(:,1) = zeros(nnod,1);u(:,2) = Sol(1:nddln:nddlt);u(:,3) = zeros(nnod,1);
+    else
+        u(:,1) = Sol(1:nddln:nddlt) ;u(:,2) = zeros(nnod,1); u(:,3) = zeros(nnod,1); 
+    end
+end
+if ndim == 2 
+    XY(:,3) = zeros(nnod,1); 
+    u(:,1) = Sol(1:nddln:nddlt);u(:,2) = Sol(2:nddln:nddlt);u(:,3) = zeros(nnod,1);
+end
+if ndim == 3  
+    u(:,1) = Sol(1:nddln:nddlt);u(:,2) = Sol(2:nddln:nddlt);u(:,3) = Sol(3:nddln:nddlt);
+end
+                 %---- calcul du facteur d'echelle et de la position deformee                    
+dX = max(XY(:,1)) - min(XY(:,1)); dY = max(XY(:,2)) - min(XY(:,2));dZ = max(XY(:,3)) - min(XY(:,3));
+                 %---- s : facteur d'echelle
+s = max([dX dY dZ])/(20.*max(abs([u(:,1) ; u(:,2) ; u(:,3)])));
+Def = XY + s * [u(:,1),u(:,2),u(:,3)];  %----- position deformee
+
+title(['deformee de la structure avec un facteur d''echelle de ',...
+        num2str(s,'%8.4f')],'Color','b')
+for iel = 1:nelt
+    loce=[];            %----- table de localisation pour l'element
+    for i=1:nnode 
+		if Connec(iel,i) > 0 loce=[loce,Connec(iel,i)]; end
+	end;  
+    Pos  = XY(loce,:);
+    X = [[Pos(:,1)]; Pos(1,1)]; Y = [[Pos(:,2)]; Pos(1,2)];Z = [[Pos(:,3)]; Pos(1,3)];
+    line(X,Y,Z,'color','g','LineStyle','--')
+    Pos = Def(loce,:);    
+    X = [[Pos(:,1)]; Pos(1,1)]; Y = [[Pos(:,2)]; Pos(1,2)];Z = [[Pos(:,3)]; Pos(1,3)];
+    line(X,Y,Z,'color','b')  
+end 
+return

Dessin/plotmodes.m

@@ -0,0 +1,36 @@
+function  plotmodes(Sol,k,f)
+% Trace des modes de vibrations d'un portique 2D  
+%
+% appel plotmodes(Sol,k,f)
+% en entree Sol : vecteur des deplacements nodaux dimension (nddlt)
+%           k   : numero du mode
+%           f   : valeur de la frequence
+%
+%  H.Oudin
+%==========================================================================
+global nddln nddlt nelt nnode
+global Coord Connec
+
+axis equal
+                 %---- calcul du facteur d'echelle et de la position deformee                    
+dX = max(Coord(:,1)) - min(Coord(:,1)); dY = max(Coord(:,2)) - min(Coord(:,2));
+%s = max([dX dY])/(10.*max(abs(Sol))) ;     %---- facteur d'echelle 
+s = max([dX dY])/(20.*max(abs([Sol(1:nddln:nddlt) ; Sol(2:nddln:nddlt)])));
+Def = Coord + s * [Sol(1:nddln:nddlt),Sol(2:nddln:nddlt)];  %----- postion deformee
+
+title([int2str(k),' mode de vibration a la frequence ',...
+        num2str(f,'%7.2f')],'Color','b')
+
+for iel = 1:nelt
+    loce=[];            %----- table de localisation pour l'element
+    for i=1:nnode 
+		if Connec(iel,i) > 0 loce=[loce,Connec(iel,i)]; end
+	end;  
+    Pos  = Coord(loce,:);
+    X = [[Pos(:,1)]; Pos(1,1)]; Y = [[Pos(:,2)]; Pos(1,2)];
+    line(X,Y,'color','g','LineStyle','--')
+    Pos = Def(loce,:);    
+    X = [[Pos(:,1)]; Pos(1,1)]; Y = [[Pos(:,2)]; Pos(1,2)];
+    line(X,Y,'color','b')   
+end 
+return

Dessin/plotstr.m

@@ -0,0 +1,66 @@
+function plotstr
+% trace du maillage d'une structure avec les numeros des noeuds et des elements 
+% et les conditions aux limites (Charges nodales, deplacements imposes) 
+%
+% appel plotstr
+%  H.Oudin
+%==========================================================================
+global nddln nnod nddlt nelt nnode ndim
+global Coord Connec Ncl F
+
+taille = get(0,'ScreenSize'); 
+figure('Name','maillage de la structure',...
+        'Position',[taille(3)/2.01 taille(4)/2.4 taille(3)/2 taille(4)/2])
+hold on
+title('maillage et CL de la structure etudiee')
+axis equal
+    grid on
+    xlabel('X');
+    ylabel('Y');
+    zlabel('Z');    
+XY = Coord;
+Cu = Ncl(1:nddln:nddlt); Cv = Ncl(2:nddln:nddlt); Cw = Ncl(3:nddln:nddlt);
+Fx = F(1:nddln:nddlt);   Fy = F(2:nddln:nddlt); Fz = F(3:nddln:nddlt);
+if ndim == 1 
+    XY(:,2) = zeros(nnod,1);XY(:,3) = zeros(nnod,1); 
+    Cv = zeros(nnod,1); Fy = zeros(nnod,1);Cw = zeros(nnod,1); Fz = zeros(nnod,1);
+end
+if ndim == 2 
+    XY(:,3) = zeros(nnod,1); 
+    Cw = zeros(nnod,1); Fz = zeros(nnod,1);
+end 
+
+
+for inod=1:nnod     % ----- Noeuds , conditions aux limites et charges nodales
+    text(XY(inod,1),XY(inod,2),XY(inod,3),[' ',int2str(inod)],'color','m','FontSize',14)
+    if Cu(inod)==1 plot3(XY(inod,1),XY(inod,2),XY(inod,3),'>','color','k','MarkerSize',10);end
+    if Cv(inod)==1 plot3(XY(inod,1),XY(inod,2),XY(inod,3),'^','color','k','MarkerSize',10);end
+    if Cw(inod)==1 plot3(XY(inod,1),XY(inod,2),XY(inod,3),'*','color','k','MarkerSize',10);end
+    if Fx(inod)~= 0  
+        text(XY(inod,1),XY(inod,2),XY(inod,3),['  \rightarrow',num2str(Fx(inod))],...
+            'HorizontalAlignment','left','color','g','FontSize',14);
+    end
+    if Fy(inod)~= 0  
+        %text(XY(inod,1),XY(inod,2),XY(inod,3),['\uparrow',num2str(Fy(inod))],...
+         %   'VerticalAlignment','bottom','color','g','FontSize',14);
+            text(XY(inod,1),XY(inod,2),XY(inod,3),['\uparrow',num2str(Fy(inod))],...
+            'VerticalAlignment','bottom','color','g','FontSize',14);
+    end   
+    if Fz(inod)~= 0  
+        text(XY(inod,1),XY(inod,2),XY(inod,3),['\uparrow',num2str(Fz(inod))],...
+            'VerticalAlignment','bottom','color','g','FontSize',14);
+    end   
+
+end
+plot3(XY(:,1),XY(:,2),XY(:,3),'.','color','m','MarkerSize',12);
+for iel = 1:nelt        %----- visualisation du maillage
+    loce=[];            %----- table de localisation pour l'element
+    for i=1:nnode 
+		if Connec(iel,i) > 0 loce=[loce,Connec(iel,i)]; end
+	end;  
+    Pos  = XY(loce,:);
+    X = [[Pos(:,1)]; Pos(1,1)];Y = [[Pos(:,2)]; Pos(1,2)];Z = [[Pos(:,3)]; Pos(1,3)];
+    line(X,Y,Z,'color','b')
+    text(mean(Pos(:,1)),mean(Pos(:,2)),mean(Pos(:,3)),[' ',int2str(iel)],'color','b','FontSize',14)
+end 
+return

Elements/Q4_ep.m

@@ -0,0 +1,60 @@
+function [Ke,Fe] = Q4_ep(iel)
+% Calcul de la matrice raideur Ke et de la force generalisee Fe 
+% pour un element Q4 d'une structure en elasticite plane
+% 
+% appel [Ke,Fe] = Q4_ep(iel)
+%    ou [Ke,Fe] = feval('Q4_ke',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire (8,8)
+%           Fe  : force generalisee elementaire (8,1)
+%
+%  H.Oudin  
+global Coord Connec Nprop Prop 
+
+%X  = Coord(Connec(iel,[1:4]),:);
+
+npg = 4;          %----- integration a 4 points de Gauss
+wg = [1,1,1,1];   %----- poids et position
+c = 1/sqrt(3); posg = [ -c -c ; c -c ; c  c ; -c  c ];
+%npg = 3; wg = [4/3, 4/3, 4/3]; posg = [sqrt(2/3) 0 ; -sqrt(1/6) sqrt(1/2) ; -sqrt(1/6) -sqrt(1/2)];
+%npg = 7; wg = [8/7, 20/63, 20/63, 20/36, 20/36, 20/36, 20/36]; 
+%c = sqrt(3/5); d = sqrt(14/5); posg = [0 0 ; 0 d ; 0 -d ; -c -c ; c -c ; c c ; -c  c ];
+
+E=Prop(Nprop(iel),1);    %----- matrice d'elasticite D
+nu=Prop(Nprop(iel),2);
+ep=Prop(Nprop(iel),3);
+if ep > 0  a = 0 ; else a = 1 ; ep = 1;  end
+coef = ep * E * (1-a*nu)/((1+nu)*(1-nu-a*nu));
+D = coef * [     1       nu/(1-a*nu)            0             ;...
+            nu/(1-a*nu)      1                  0             ;...
+                 0           0       .5*(1-nu-a*nu)/(1-a*nu)];
+
+ndle = 8;            %----- initialisations
+Ke = zeros(ndle); Fe = zeros(ndle,1); %  aire=0
+for ipg=1:npg       %----- boucle d'integration
+   s = posg(ipg,1); t = posg(ipg,2); poids = wg(ipg);
+    %----- vecteur <N(s,t)>
+   N = .25*[(1-s)*(1-t)  (1+s)*(1-t)  (1+s)*(1+t)  (1-s)*(1+t)]; 
+    %----- matrice [dN/ds ;dN/dt]
+   dN = .25*[-(1-t)  (1-t) (1+t)  -(1+t)
+   		     -(1-s) -(1+s) (1+s)   (1-s)]; 
+	%----- matrice jacobienne
+   J = dN*Coord(Connec(iel,[1:4]),:);
+   detj = J(1,1)*J(2,2)-J(1,2)*J(2,1);
+   J_1 = [J(2,2) -J(1,2); -J(2,1) J(1,1)]/detj ;
+    %----- matrice [dN/dx ;dN/dy]
+   dNx = J_1*dN;
+	%----- matrice B(3x8)
+   B=zeros(3,8);
+   B(1,[1 3 5 7])=dNx(1,:);
+   B(2,[2 4 6 8])=dNx(2,:);			
+   B(3,[1 3 5 7,2 4 6 8])=[dNx(2,:),dNx(1,:)];
+	%----- matrice Ke(8x8)   %aire=aire+detj*poids;  
+   Ke=Ke+(B'*D*B)*detj*poids;
+	%----- vecteur Fe(8,1) 
+   fx=Prop(Nprop(iel),4); fy=Prop(Nprop(iel),5);
+   Fe([1 3 5 7],1) = Fe([1 3 5 7],1)+ ep*fx*detj*poids*N';
+   Fe([2 4 6 8],1) = Fe([2 4 6 8],1)+ ep*fy*detj*poids*N';
+end
+%disp(Ke),disp(Fe)
+return

Elements/Q4_ep_stress.m

@@ -0,0 +1,59 @@
+function sig = Q4_ep_stress(iel,Ue)
+% Calcul du tenseur "contrainte moyenne" sur un element Q4
+% en elasticite plane  
+% 
+% appel [sig] = Q4_ep_stress(iel,Ue)
+%    ou [sig] = feval('Q4_stress',iel,Ue)
+% en entree iel : numero de l'element
+%           Ue  : vecteur des depalements nodaux de l'element
+% en sortie Tenseur des contraintes moyennes sur l'element
+%
+%  H.Oudin 
+global Coord Connec Nprop Prop 
+
+                  %-----  position des points de calcul de la contrainte
+npg = 4; c = 1/sqrt(3); posg = [ -c -c ; c -c ; c  c ; -c  c ]; %points de gauss
+%npg = 4;  c = 1;         posg = [ -c -c ; c -c ; c  c ; -c  c ]; %noeuds de l'element
+%npg = 1; posg = [ 0 0 ]; %au centre
+
+
+E=Prop(Nprop(iel),1);    %----- matrice d'elasticite D
+nu=Prop(Nprop(iel),2);
+ep=Prop(Nprop(iel),3);
+if ep > 0  a = 0 ; else a = 1 ; ep = 1;  end
+coef = E * (1-a*nu)/((1+nu)*(1-nu-a*nu));
+D = coef * [     1       nu/(1-a*nu)            0             ;...
+            nu/(1-a*nu)      1                  0             ;...
+                 0           0       .5*(1-nu-a*nu)/(1-a*nu)];
+
+ndle = 8;            %----- initialisations
+sigpg = zeros(4,3);sig = zeros(1,3);
+
+%fprintf('Valeur des contraintes en des points de l''element %3i\n',iel)
+for ipg=1:npg       %----- boucle sur les points 
+   s = posg(ipg,1); t = posg(ipg,2);
+    %----- vecteur <N(s,t)>
+   N = .25*[(1-s)*(1-t)  (1+s)*(1-t)  (1+s)*(1+t)  (1-s)*(1+t)]; 
+    %----- matrice [dN/ds ;dN/dt]
+   dN = .25*[-(1-t)  (1-t) (1+t)  -(1+t)
+   		     -(1-s) -(1+s) (1+s)   (1-s)]; 
+	%----- matrice jacobienne
+   J = dN*Coord(Connec(iel,[1:4]),:);
+   detj = J(1,1)*J(2,2)-J(1,2)*J(2,1);
+   J_1 = [J(2,2) -J(1,2); -J(2,1) J(1,1)]/detj ;
+    %----- matrice [dN/dx ;dN/dy]
+   dNx = J_1*dN;
+	%----- matrice B(3x8)
+   B=zeros(3,8);
+   B(1,[1 3 5 7])=dNx(1,:);
+   B(2,[2 4 6 8])=dNx(2,:);			
+   B(3,[1 3 5 7,2 4 6 8])=[dNx(2,:),dNx(1,:)];
+	%----- tableau des contraintes aux points de Gauss   
+   sigpg(ipg,:) = (D*B*Ue)';
+   sig(1,:) = sig(1,:) + sigpg(ipg,:);
+  % fprintf('point (s,t) %5.2f %5.2f ',s,t)
+  % fprintf('sigxx = %8.3e  sigyy = %8.3e   sigxy = %8.3e  \n',sigpg(ipg,1),sigpg(ipg,2),sigpg(ipg,3))
+end
+%----- Valeur moyenne du tenseur des contraintes
+sig = sig/npg;
+return

Elements/Q4_th.m

@@ -0,0 +1,46 @@
+function [Ke,Fe] = Q4_th(iel)
+% Calcul de la matrice comportement Ke et du vecteur flux Fe 
+% pour un element Q4 de thermique lineaire en regime permanent
+% 
+% appel [Ke,Fe] = Q4_ke(iel)
+%    ou [Ke,Fe] = feval('Q4_ke',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire (4,4)
+%           Fe  : force generalisee elementaire (4,1)
+%
+%  H.Oudin  
+global Coord Connec Nprop Prop 
+
+npg = 4;          %----- integration a 4 points de Gauss
+wg = [1,1,1,1];   %----- poids et position
+c = 1/sqrt(3); posg = [ -c -c ; c -c ; c  c ; -c  c ];
+
+D=Prop(Nprop(iel),1);    %----- carateristiques thermiques
+r=Prop(Nprop(iel),2);
+
+ndle = 4;            %----- initialisations
+Ke = zeros(ndle); Fe = zeros(ndle,1); %  aire=0
+for ipg=1:npg       %----- boucle d'integration
+   s = posg(ipg,1); t = posg(ipg,2); poids = wg(ipg);
+    %----- vecteur <N(s,t)>
+   N = .25*[(1-s)*(1-t)  (1+s)*(1-t)  (1+s)*(1+t)  (1-s)*(1+t)]; 
+    %----- matrice [dN/ds ;dN/dt]
+   dN = .25*[-(1-t)  (1-t) (1+t)  -(1+t)
+   		     -(1-s) -(1+s) (1+s)   (1-s)]; 
+	%----- matrice jacobienne
+   J = dN*Coord(Connec(iel,[1:4]),:);
+   detj = J(1,1)*J(2,2)-J(1,2)*J(2,1);
+   J_1 = [J(2,2) -J(1,2); -J(2,1) J(1,1)]/detj ;
+    %----- matrice [dN/dx ;dN/dy]
+   dNx = J_1*dN;
+	%----- matrice B(2x4)
+   B=zeros(2,4);
+   B(1,[1 2 3 4])=dNx(1,:);
+   B(2,[1 2 3 4])=dNx(2,:);			
+	%----- matrice Ke(4x4)   %aire=aire+detj*poids;  
+   Ke=Ke+(B'*D*B)*detj*poids;
+	%----- vecteur Fe(4,1) 
+   Fe([1 2 3 4],1) = Fe([1 2 3 4],1)+ r*detj*poids*N';
+end
+%disp(Ke),disp(Fe)
+return

Elements/T3_ep.m

@@ -0,0 +1,47 @@
+function [Ke,Fe] = T3_ep(iel)
+% Calcul de la matrice raideur Ke et de la force generalisee Fe 
+% pour un element T3 d'une structure en elasticite plane
+% 
+%   Nous n'utilisons pas l'integration numerique, le calcul des
+%   matrices est presente dans le cours chapitre V-2.4 nous utilisons
+%   ici ces resultats.
+%
+% appel [Ke,Fe] = T3_ke(iel)
+%    ou [Ke,Fe] = feval('T3_ke',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire (6,6)
+%           Fe  : force generalisee elementaire (6,1)
+%
+%  H.Oudin 
+global Coord Connec Nprop Prop 
+
+ndle = 6;                %----- initialisations
+Ke = zeros(ndle); Fe = zeros(ndle,1); 
+
+E=Prop(Nprop(iel),1);    %----- matrice d'elasticite D
+nu=Prop(Nprop(iel),2);
+ep=Prop(Nprop(iel),3);
+if ep > 0  a = 0 ; else a = 1 ; ep = 1;  end
+coef = ep * E * (1-a*nu)/((1+nu)*(1-nu-a*nu));
+D = coef * [     1       nu/(1-a*nu)            0             ;...
+            nu/(1-a*nu)      1                  0             ;...
+                 0           0       .5*(1-nu-a*nu)/(1-a*nu)];
+
+X  = Coord(Connec(iel,[1:3]),:);    % coordonnees des noeuds de l'element
+
+                                    %-----  determinant detj = 2A
+detj=(X(2,1)-X(1,1))*(X(3,2)-X(1,2))-(X(3,1)-X(1,1))*(X(2,2)-X(1,2));
+        
+B=zeros(3,6);                       %-----  matrice B(3,6) 
+B(1,[1 3 5])=(X([2 3 1],2)-X([3 1 2],2) )' /detj;
+B(2,[2 4 6])=(X([3 1 2],1)-X([2 3 1],1) )' /detj ;
+B(3,[1 3 5 2 4 6]) =[ B(2,[2 4 6]), B(1,[1 3 5]) ];
+        
+Ke = (B'*D*B)*.5*detj;              %----- matrice Ke(6,6)
+                                    %----- vecteur Fe(6,1)
+fx=Prop(Nprop(iel),4); fy=Prop(Nprop(iel),5);
+Fe([1 3 5],1) = (ep*fx*detj/6)*[1 1 1]';
+Fe([2 4 6],1) = (ep*fy*detj/6)*[1 1 1]';
+
+%disp(Ke),disp(Fe)
+return

Elements/T3_ep_stress.m

@@ -0,0 +1,36 @@
+function sig = T3_ep_stress(iel,Ue)
+% Calcul du champ de contrainte sur un element T3
+% en elasticite plane  (c'est une constante)
+% 
+% appel [sig] = T3_ep_stress(iel,Ue)
+%    ou [sig] = feval('T3_stress',iel,Ue)
+% en entree iel : numero de l'element
+%           Ue  : vecteur des depalements nodaux de l'element
+% en sortie Tenseur des contraintes moyennes sur l'element
+%
+%  H.Oudin 
+global Coord Connec Nprop Prop 
+
+nnode=3; ndle = 6;                %----- initialisations
+
+E=Prop(Nprop(iel),1);    %----- matrice d'elasticite D
+nu=Prop(Nprop(iel),2);
+ep=Prop(Nprop(iel),3);
+if ep > 0  a = 0 ; else a = 1 ; ep = 1;  end
+coef = E * (1-a*nu)/((1+nu)*(1-nu-a*nu));
+D = coef * [     1       nu/(1-a*nu)            0             ;...
+            nu/(1-a*nu)      1                  0             ;...
+                 0           0       .5*(1-nu-a*nu)/(1-a*nu)];
+
+X  = Coord(Connec(iel,[1:nnode]),:);    % coordonnees des noeuds de l'element
+
+                                    %-----  determinant detj = 2A
+detj=(X(2,1)-X(1,1))*(X(3,2)-X(1,2))-(X(3,1)-X(1,1))*(X(2,2)-X(1,2));
+        
+B=zeros(nnode,ndle);                       %-----  matrice B(3,6) 
+B(1,[1 3 5])=(X([2 3 1],2)-X([3 1 2],2) )' /detj;
+B(2,[2 4 6])=(X([3 1 2],1)-X([2 3 1],1) )' /detj ;
+B(3,[1 3 5 2 4 6]) =[ B(2,[2 4 6]), B(1,[1 3 5]) ];
+        
+sig = (D*B*Ue)';              %----- vecteur des contraintes
+return

Elements/barre_ke.m

@@ -0,0 +1,40 @@
+function [Ke,Fe] = barre_ke(iel)
+% Calcul de la matrice raideur Ke et de la force generalisee Fe 
+% pour un element (iel) d'une structure treillis
+%
+% appel [Ke,Fe] = barre_ke(iel)
+%    ou [Ke,Fe] = feval('barre_ke',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire (2*ndim,2*ndim)
+%           Fe  : force generalisee elementaire (2*ndim) 
+%
+%  H.Oudin 
+global ndim
+global Coord Connec Nprop Prop 
+
+ES=Prop(Nprop(iel),1);
+X  = Coord(Connec(iel,:),:);    % ----- Coordonnees des 2 noeuds de l'element
+dX = X(2,:) - X(1,:);           % ----- x2-x1  et y2-y1
+if ndim == 1
+    L = abs(dX);
+    Ke = (ES/L)*[  1  -1
+                  -1   1 ];
+    f=Prop(Nprop(iel),2);
+    Fe = (f*L/2)*[1;1];
+elseif ndim == 2
+    L = sqrt(dX(1)^2 + dX(2)^2);
+    c = dX(1)/L; s = dX(2)/L;   % ----- Cosinus directeurs de l'element
+    cc = c*c; ss = s*s; cs = c*s;
+    Ke = (ES/L)*[  cc  cs -cc -cs
+                   cs  ss -cs -ss
+                  -cc -cs  cc  cs
+                  -cs -ss  cs  ss];
+     fx=Prop(Nprop(iel),2); fy=Prop(Nprop(iel),3);
+     Fe = (L/2)*[fx;fy;fx;fy];
+     %Fe =(L/2)*[fx*c-fy*s; fx*s+fy*c; fx*c-fy*s; fx*s+fy*c];
+ elseif ndim == 3
+     disp('================================================ ');
+     disp('       element non programme ');
+     disp('================================================ ');
+ end
+return

Elements/barre_keme.m

@@ -0,0 +1,43 @@
+function [Ke,Me] = barre_keme(iel)
+% Calcul de la matrice raideur Ke et de la matrice masse 
+% pour un element (iel) d'une structure treillis
+%
+% appel [Ke,Me] = barre_keme(iel)
+%    ou [Ke,Me] = feval('barre_keme',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire (2*ndim,2*ndim)
+%           Me  : matrice masse   elementaire (2*ndim,2*ndim)
+%
+%  H.Oudin 
+global ndim
+global Coord Connec Nprop Prop 
+
+ES=Prop(Nprop(iel),1);
+roS=Prop(Nprop(iel),2);
+X  = Coord(Connec(iel,:),:);    % ----- Coordonnees des 2 noeuds de l'element
+dX = X(2,:) - X(1,:);           % ----- x2-x1  et y2-y1
+if ndim == 1
+    L = abs(dX);
+    Ke = (ES/L)*[  1  -1
+                  -1   1 ];  
+    Me = (roS*L/6)*[2  1
+                    1  2 ];
+elseif ndim == 2
+    L = sqrt(dX(1)^2 + dX(2)^2);
+    c = dX(1)/L; s = dX(2)/L;   % ----- Cosinus directeurs de l'element
+    cc = c*c; ss = s*s; cs = c*s;
+    Ke = (ES/L)*[  cc  cs -cc -cs
+                   cs  ss -cs -ss
+                  -cc -cs  cc  cs
+                  -cs -ss  cs  ss];
+    M = (roS*L/6)*[2  1
+                   1  2 ];
+    P = [  c  s  0  0  
+            0  0  c  s ];
+    Me =P' * M * P;            
+ elseif ndim == 3
+     disp('================================================ ');
+     disp('       element non programme ');
+     disp('================================================ ');
+ end
+return

Elements/barre_stress.m

@@ -0,0 +1,28 @@
+function Ne = barre_stress(iel,Ue)
+% Calcul de l'effort normal dans un element d'une structure treillis
+%
+% appel Ne=feval('barre_stress',iel,Ue)
+% en entree iel : numero de l'element
+%           Ue  : vecteur des depalements nodaux de l'element
+% en sortie l'effort normal dans la barre
+% H.Oudin 
+global ndim
+global Coord Connec Nprop Prop 
+
+ES=Prop(Nprop(iel),1);
+X  = Coord(Connec(iel,:),:);
+dX = X(2,:) - X(1,:);
+
+if ndim == 1
+    L = abs(dX);
+    Ne = (ES/L)*(Ue(2)-Ue(1));
+elseif ndim == 2
+    L = sqrt(dX(1)^2 + dX(2)^2);
+    c = dX(1)/L; s = dX(2)/L; 
+    Ne=(ES/L)*(c*(Ue(3)-Ue(1))+s*(Ue(4)-Ue(2)));
+ elseif ndim == 3
+     disp('================================================ ');
+     disp('       element non programme ');
+     disp('================================================ ');
+end
+return

Elements/poutre_ke.m

@@ -0,0 +1,107 @@
+function [Ke,Fe] = poutre_ke(iel)
+% Calcul de la matrice raideur Ke et de la force generalisee Fe 
+% pour un element (iel) d'une structure portique
+% 
+% appel [Ke,Fe] = poutre_ke(iel)
+%    ou [Ke,Fe] = feval('poutre_ke',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire 
+%           Fe  : force generalisee elementaire  
+%
+%  H.Oudin & D.Motte (EF 3D)
+global ndim
+global Coord Connec Nprop Prop 
+
+X  = Coord(Connec(iel,:),:);
+dX = X(2,:) - X(1,:);  
+
+if ndim == 1  
+    EI=Prop(Nprop(iel),1);
+    L = abs(dX);    %----- matrice raideur elementaire sur (v,teta)
+    Ke = (EI/L^3)*[  12   6*L   -12   6*L
+                     6*L  4*L^2 -6*L  2*L^2
+                    -12  -6*L    12  -6*L
+                     6*L  2*L^2 -6*L  4*L^2
+                   ]; 
+                   
+    f=Prop(Nprop(iel),2);
+    Fe = (f*L/2)*[1; L/6; 1; -L/6];
+elseif ndim == 2 
+ ES=Prop(Nprop(iel),1); EI=Prop(Nprop(iel),2);
+ L = sqrt(dX(1)^2 + dX(2)^2);
+ c = dX(1)/L; s = dX(2)/L;  
+ K = (ES/L)*[  1  0  0 -1  0  0
+               0  0  0  0  0  0
+               0  0  0  0  0  0
+              -1  0  0  1  0  0
+               0  0  0  0  0  0
+               0  0  0  0  0  0
+            ];
+ K = K + (EI/L^3)*[  0  0    0     0   0    0
+                     0  12   6*L   0  -12   6*L
+                     0  6*L  4*L^2 0  -6*L  2*L^2
+                     0  0    0     0   0    0
+                     0 -12  -6*L   0   12  -6*L
+                     0  6*L  2*L^2 0  -6*L  4*L^2
+                   ]; 
+  P = [  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
+            ];
+  Ke =P' * K * P;    %----- matrice raideur elementaire sur (u,v,teta)
+
+  f  = [ c*Prop(Nprop(iel),3)+s*Prop(Nprop(iel),4)
+        -s*Prop(Nprop(iel),3)+c*Prop(Nprop(iel),4)];
+  fe = (f(1)*L/2)*[1; 0; 0; 1; 0; 0];
+  fe = fe + (f(2)*L/2)*[0; 1; L/6; 0; 1; -L/6];
+  Fe = P' * fe;     %----- force generalisee elementaire sur (u,v,teta)
+  
+elseif ndim == 3
+%-------------------------------------------------------------------------------------------------
+ES = Prop(Nprop(iel),1);				% Matrice de raideur locale
+GJ = Prop(Nprop(iel),2);
+EIy = Prop(Nprop(iel),3);
+EIz = Prop(Nprop(iel),4);
+L = sqrt(dX(1)^2+dX(2)^2+dX(3)^2);
+
+K = [ ES/L	0			0			0		0			0			-ES/L	0			0			0		0			0
+      0		12*EIz/L^3	0			0		0			6*EIz/L^2	0		-12*EIz/L^3	0			0		0			6*EIz/L^2
+	  0		0			12*EIy/L^3	0		-6*EIy/L^2	0			0		0			-12*EIy/L^3	0		-6*EIy/L^2	0
+	  0		0			0			GJ/L	0			0			0		0			0			-GJ/L	0			0
+	  0		0			-6*EIy/L^2	0		4*EIy/L		0			0		0			6*EIy/L^2	0		2*EIy/L		0
+      0		6*EIz/L^2	0			0		0			4*EIz/L		0		-6*EIz/L^2	0			0		0			2*EIz/L
+	  -ES/L	0			0			0		0			0			ES/L	0			0			0		0			0
+      0		-12*EIz/L^3	0			0		0			-6*EIz/L^2	0		12*EIz/L^3	0			0		0			-6*EIz/L^2
+	  0		0			-12*EIy/L^3	0		6*EIy/L^2	0			0		0			12*EIy/L^3	0		6*EIy/L^2	0
+	  0		0			0			-GJ/L	0			0			0		0			0			GJ/L	0			0
+	  0		0			-6*EIy/L^2	0		2*EIy/L		0			0		0			6*EIy/L^2	0		4*EIy/L		0
+      0		6*EIz/L^2	0			0		0			2*EIz/L		0		-6*EIz/L^2	0			0		0			4*EIz/L
+	];
+
+tx=dX(1)/L;								% parametres du repere local / repere global
+ty=dX(2)/L;
+tz=dX(3)/L;
+t=[tx ; ty ; tz];
+prodzt=cross([0,0,1],t);
+tn=prodzt/sqrt(prodzt(1)^2+prodzt(2)^2+prodzt(3)^2);
+tm=cross(t,tn);
+Q=[ t(1)	,	t(2)	,	t(3)		% Matrice de passage 3x3
+	tn(1)	,	tn(2)	,	tn(3)
+	tm(1)	,	tm(2)	,	tm(3)
+  ];
+nul=zeros(3);
+P=[	Q	,	nul	,	nul	,	nul			% Matrice de passage 12x12
+	nul	,	Q	,	nul	,	nul
+	nul	,	nul	,	Q	,	nul
+	nul	,	nul	,	nul	,	Q
+  ];
+
+Ke = P' * K * P;		%----- matrice raideur elementaire sur (u,v,w,teta1,teta2,teta3)
+
+Fe = zeros(12,1);		%----- pas de force generalisee elementaire
+%-------------------------------------------------------------------------------------------------
+end
+return

Elements/poutre_keme.m

@@ -0,0 +1,72 @@
+function [Ke,Me] = poutre_keme(iel)
+% Calcul des matrices elementaire Ke et Me de l'element (iel) d'un portique
+% 
+% appel [Ke,Me] = poutre_keme(iel)
+%    ou [Ke,Me] = feval('poutre_keme',iel)
+% en entree iel : numero de l'element
+% en sortie Ke  : matrice raideur elementaire (4,4)en 1D; (6,6) en 2D 
+%           Me  : matrice masse elementaire 
+%
+%  H.Oudin 
+global ndim
+global Coord Connec Nprop Prop 
+
+X  = Coord(Connec(iel,:),:);
+dX = X(2,:) - X(1,:);  
+
+if ndim == 1  
+    EI=Prop(Nprop(iel),2);RS=Prop(Nprop(iel),3);
+    L = abs(dX);    %----- matrice raideur elementaire sur (v,teta)
+    Ke = (EI/L^3)*[  12   6*L   -12   6*L
+                     6*L  4*L^2 -6*L  2*L^2
+                    -12  -6*L    12  -6*L
+                     6*L  2*L^2 -6*L  4*L^2
+                   ]; 
+                   
+     Me = (RS*L)*[  13/35     11*L/210    9/70   -13*L/420
+                    11*L/210   L*L/105  13*L/420  -L*L/140 
+                    9/70      13*L/420    13/35   -11*L/210
+                  -13*L/420   -L*L/140  -11*L/210  L*L/105
+               ];
+elseif ndim == 2 
+ ES=Prop(Nprop(iel),1); EI=Prop(Nprop(iel),2);RS=Prop(Nprop(iel),3);
+ L = sqrt(dX(1)^2 + dX(2)^2);
+ c = dX(1)/L; s = dX(2)/L;  
+ K = (ES/L)*[  1  0  0 -1  0  0
+               0  0  0  0  0  0
+               0  0  0  0  0  0
+              -1  0  0  1  0  0
+               0  0  0  0  0  0
+               0  0  0  0  0  0
+            ];
+ K = K + (EI/L^3)*[  0  0    0     0   0    0
+                     0  12   6*L   0  -12   6*L
+                     0  6*L  4*L^2 0  -6*L  2*L^2
+                     0  0    0     0   0    0
+                     0 -12  -6*L   0   12  -6*L
+                     0  6*L  2*L^2 0  -6*L  4*L^2
+                   ]; 
+  P = [  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
+            ];
+  Ke =P' * K * P;    %----- matrice raideur elementaire sur (u,v,teta)
+
+  M = (RS*L)*[  1/3      0            0      1/6     0         0
+                 0     13/35     11*L/210     0     9/70    -13*L/420
+                 0    11*L/210      L*L/105   0    13*L/420  -L*L/140 
+                1/6      0             0     1/3     0          0
+                 0      9/70      13*L/420    0     13/35    -11*L/210
+                 0   -13*L/420     -L*L/140   0    -11*L/210   L*L/105
+               ];
+  Me =P' * M * P;    %----- matrice massse elementaire sur (u,v,teta)
+  
+elseif ndim == 3
+     disp('================================================ ');
+     disp('       element non programme ');
+     disp('================================================ ');
+end
+return

Elements/poutre_stress.m

@@ -0,0 +1,46 @@
+function [Ne,Te,Mfe] = poutre_stress(iel,Ue)
+% Calcul des contraintes dans un element poutre d'une structure portique
+%
+% appel poutre_stress(iel,Ue)
+%    ou feval('poutre_stress',iel,Ue)
+% en entree iel : numero de l'element
+%           Ue  : vecteur des depalements nodaux de l'element
+%
+%  H.Oudin 
+global ndim
+global Coord Connec Nprop Prop 
+
+X  = Coord(Connec(iel,:),:);
+dX = X(2,:) - X(1,:);
+
+if ndim == 1
+    EI=Prop(Nprop(iel),1);
+    L = abs(dX); Ne = 0;
+    Te = -(EI/L^3)*(12*Ue(1)+6*L*Ue(2)-12*Ue(3)+6*L*Ue(4));
+    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)];
+    fprintf('Dans l''element %3i\n',iel)
+    fprintf('Ne = %8.3e  Te = %8.3e Mfe au noeud i  %8.3e  en noeud j %8.3e \n',Ne,Te,Mfe(1),Mfe(2))
+elseif ndim == 2
+    ES=Prop(Nprop(iel),1); EI=Prop(Nprop(iel),2);
+    L = sqrt(dX(1)^2 + dX(2)^2);
+    c = dX(1)/L; s = dX(2)/L; 
+    P = [  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
+         ];
+    Ue = P*Ue;
+    Ne = (ES/L)*(Ue(4)-Ue(1));
+    Te = -(EI/L^3)*(12*Ue(2)+6*L*Ue(3)-12*Ue(5)+6*L*Ue(6));
+    Mfe= (EI/L^2)*[-6*Ue(2)-4*L*Ue(3)+6*Ue(5)-2*L*Ue(6) ; 6*Ue(2)+2*L*Ue(3)-6*Ue(5)+4*L*Ue(6)];
+    fprintf('Dans l''element %3i\n',iel)
+    fprintf('Ne = %8.3e  Te = %8.3e Mfe au noeud i  %8.3e  en noeud j %8.3e \n',Ne,Te,Mfe(1),Mfe(2))
+ elseif ndim == 3
+     disp('================================================ ');
+     disp('       element non programme ');
+     disp('================================================ ');
+end
+
+return

Generaux/resultante.m

@@ -0,0 +1,23 @@
+function [Fx,Fy,Fz] = resultante(F)
+% Calcul des resultantes d'un vecteur F 
+% en fonction de la dimension 1D 2D ou 3D du probleme
+%
+% appel [Fx,Fy,Fz] = resultante(F)
+%    ou [Fx,Fy,Fz] = feval('resultante',F)
+% en entree F   : vecteur de dimension nddlt
+% en sortie les resltantes Fx, 0, O en dim 1
+%                          Fx,Fy, 0 en dim 2
+%                          Fx,Fy,Fz en dim 3
+% 
+%  H.Oudin 
+global nddln nddlt ndim
+
+if ndim == 3     
+    Fz=sum(F(3:nddln:nddlt)); Fy=sum(F(2:nddln:nddlt)); Fx=sum(F(1:nddln:nddlt));
+elseif ndim == 2 
+    Fz = 0 ; Fy=sum(F(2:nddln:nddlt)); Fx=sum(F(1:nddln:nddlt));
+else             
+    Fz = 0 ; Fy = 0 ; Fx=sum(F(1:nddln:nddlt));
+end
+
+return

Sol_analytique/barre_compar.m

@@ -0,0 +1,36 @@
+function  barre_compar(U)
+% comparaison de la solution numerique et analytique pour 
+% une colonne soumise à son poids propre
+%  H.Oudin 
+global nddln nelt nnode
+global Coord Connec Nprop Prop
+
+fa1=@(x) -(6*x/1000).*(1-x/1200);  % solution analytique
+fa2=@(x) -6*(1-x/600);
+
+taille = get(0,'ScreenSize'); 
+figure('Name','comparaison avec la solution analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.4 taille(3)/2 taille(4)/2])          
+           %----- solution analytique pour l'exemple du chapitre II_3 du cours
+subplot(2,1,1), title('deplacement u'), hold on,
+fplot(fa1,[0 600],'r')
+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);
+    X = Coord(Connec(iel,:)); x2 = X(2); x1=X(1); L = abs(x2-x1);  
+    x=x1:(x2-x1)/5: x2; y=Ue(1)+(Ue(2)-Ue(1))*(x-x1)/L; 
+    plot(x,y,'b')
+end
+grid
+subplot(2,1,2), title('effort normal N'), hold on 
+fplot(fa2,[0 600],'r') , 
+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);
+    X = Coord(Connec(iel,:)); x2 = X(2); x1=X(1); L = abs(x2-x1);  
+    ES=Prop(Nprop(iel),1); Ne=(ES/L)*(Ue(2)-Ue(1));
+    x=[x1:(x2-x1)/5:x2];
+    y=Ne*ones(1,length(x));
+    plot(x,y,'b')
+end
+return

Sol_analytique/poutre_compar.m

@@ -0,0 +1,42 @@
+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

Sol_analytique/vibra_barre_exo2.m

@@ -0,0 +1,68 @@
+clc;clear all;close all;
+% Solution analytique frequence et mode de vibration d'une barre 
+% Barre bi-encastree ;libre -libre ; encastree - ressort
+%  H.Oudin 
+%----------------------------------------------------------
+% Caracteristiques de la barre
+ES= 1 ; roS=1; L= 1;
+%----------------------------------------------------------
+k = 0;
+taille = get(0,'ScreenSize'); 
+
+while k < 5
+k = menu('Choix des conditions aux limites',...
+    'bi-encastree','encastree - libre','libre -libre',...
+    'encastree - ressort', 'sortir') ;
+tr = [1,2,3,4,5];
+
+nmod = 3;
+ switch tr(k)
+  case 1, disp('Barre bi-encastree'); 
+  figure('Name','solution analytique : 3 premiers modes de vibration d''une barre',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+    for imod = 1:nmod
+     anal=imod*pi*sqrt(ES/roS/L/L);
+     subplot(nmod,1,imod), fplot(@(x) sin(imod*pi*x/L),[0 L],'r'),...
+     title(['barre bi-encastree : ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+    end     
+  case 2, disp('encastree -libre');
+  figure('Name','solution analytique : 3 premiers modes de vibration d''une barre',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+    for imod = 1:nmod
+     aa=(2*imod-1)*pi/2;
+     anal=aa*sqrt(ES/roS/L/L);
+     subplot(nmod,1,imod), fplot(@(x) sin(aa*x/L),[0 L],'r'),...
+     title(['barre encastree -libre : ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+    end     
+  case 3, disp('libre - libre');
+  figure('Name','solution analytique : 3 premiers modes de vibration d''une barre',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+    for imod = 1:nmod
+     anal=imod*pi*sqrt(ES/roS/L/L);
+     subplot(nmod,1,imod), fplot(@(x) cos(imod*pi*x/L),[0 L],'r'),...
+     title(['barre libre -libre : ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+    end
+  case 4, disp('encastree - ressort');
+  alpa = input('donner la valeur du rapport ES/kl ? [1/2]: ');
+  if isempty(alpa) alpa=0.5; end 
+  figure('Name',['equation implicite tan(x)+alpa*x = 0  pour alpa =',...
+      num2str(alpa)],'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+  subplot(nmod+1,1,1),hold on, fplot(@(x) tan(x),[0,3*pi,-6,6],'r'),...
+           fplot(@(x) -alpa*x,[0,3*pi,-6,6],'g'),
+       title('representation de l''equation implicite'), grid
+ %figure('Name','solution analytique : 3 premiers modes de vibration d''une barre',...
+%        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])      
+  for imod = 1:nmod
+     a= (2*imod-1)*pi/2+0.01;
+     b=(2*imod+1)*pi/2-0.01;
+     lx = fzero(@(x) tan(x)+alpa*x,[a b]);
+     anal=lx*sqrt(ES/roS/L/L);
+     subplot(nmod+1,1,imod+1), fplot(@(x) sin(lx*x/L),[0 L],'r'),...
+     title(['barre encastree - ressort : ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+    end
+ end
+end

Sol_analytique/vibra_poutre.m

@@ -0,0 +1,100 @@
+clc;clear all;close all;
+% Solution analytique frequence et mode de vibration d'une poutre 
+%  H.Oudin 
+%----------------------------------------------------------
+% Caracteristiques de la poutre
+EI= 1 ; roS=1; L= 1;
+%----------------------------------------------------------
+k = 0;
+taille = get(0,'ScreenSize'); 
+
+while k < 6
+k = menu('Choix des conditions aux limites',...
+    'libre -libre','bi-appuyee','bi-encastree',...
+    'appuyee - libre','encastree - appuyee', 'sortir') ;
+tr = [1,2,3,4,5,6];
+nmod = 3;
+switch tr(k)
+case 1, disp('poutre libre - libre');
+  disp('c''est l''exemple traite dans le cours');
+  disp('il existe deux modes rigide une translation une rotation');
+  figure('Name','3 premiers modes de vibration d''une poutre libre - libre',...
+  'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+  subplot(nmod+1,1,1),hold on, fplot('[cos(x) , 1/cosh(x)]',[0,5*pi,-1,1],'b'),...
+  title('representation de l''equation implicite cos(x)cosh(x)-1 = 0'), grid   
+  for imod = 1:nmod
+     a= imod*pi;
+     b=(imod+1)*pi;
+     lx = fzero(@(x) cos(x)*cosh(x)-1,[a b])
+     anal=(lx^2)*sqrt(EI/roS/(L^4));
+     aa = -(cos(lx)-cosh(lx))/(sin(lx)-sinh(lx));
+     mode = @(x) cos(lx*x/L)+cosh(lx*x/L)+aa*(sin(lx*x/L)+sinh(lx*x/L));
+     %mode2= @(x) (cos(lx*x/L)+cosh(lx*x/L)+aa*(sin(lx*x/L)+sinh(lx*x/L))).*...
+     %(cos(lx*x/L)+cosh(lx*x/L)+aa*(sin(lx*x/L)+sinh(lx*x/L)));
+     %m = roS*quad(mode2,0,L);  % Les modes sont M norme  m=1 
+     subplot(nmod+1,1,imod+1), fplot(mode,[0 L],'r'),...
+     title(['poutre libre - libre ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+   end  
+case 2, disp('poutre bi-appuyee les solutions sont lx = i*pi'); 
+    figure('Name','3 premiers modes de vibration d''une poutre bi-appuyee',...
+       'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+    for imod = 1:nmod
+     lx = imod*pi
+     anal=(lx^2)*sqrt(EI/roS/(L^4));
+     subplot(nmod,1,imod), fplot(@(x) sin(imod*pi*x/L),[0 L],'r'),...
+     title(['poutre bi-appuyee ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+    end     
+case 3, disp('poutre bi-encastree');
+  figure('Name','3 premiers modes de vibration d''une poutre bi-encastree',...
+  'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+  subplot(nmod+1,1,1),hold on, fplot('[cos(x) , 1/cosh(x)]',[0,5*pi,-1,1],'b'),...
+  title('representation de l''equation implicite cos(x)cosh(x)-1 = 0'), grid    
+  for imod = 1:nmod
+     a= imod*pi;
+     b=(imod+1)*pi;
+     lx = fzero(@(x) cos(x)*cosh(x)-1,[a b])
+     anal=(lx^2)*sqrt(EI/roS/(L^4));
+     aa = -(cos(lx)-cosh(lx))/(sin(lx)-sinh(lx));
+     mode = @(x) cos(lx*x/L)-cosh(lx*x/L)+aa*(sin(lx*x/L)-sinh(lx*x/L));
+     subplot(nmod+1,1,imod+1), fplot(mode,[0 L],'r'),...
+     title(['poutre bi - encastree ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+  end  
+case 4, disp('poutre appuyee - libre');
+  disp('il existe un mode rigide de rotation');
+  figure('Name','3 premiers modes de vibration d''une poutre appuyee - libre',...
+  'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+  subplot(nmod+1,1,1),hold on, fplot('[tan(x) , tanh(x)]',[0,4*pi,-1.5,1.5],'b'),...
+  title('representation de l''equation implicite tan(x) = th(x)'), grid   
+  for imod = 1:nmod
+     a= imod*pi/2+.1;
+     b=(2*imod+1)*pi/2-.1;
+     lx = fzero(@(x) tan(x)-tanh(x),[a b])
+     anal=(lx^2)*sqrt(EI/roS/(L^4));
+     aa = sin(lx)/sinh(lx);
+     mode = @(x) sin(lx*x/L)+aa*sinh(lx*x/L);
+     subplot(nmod+1,1,imod+1), fplot(mode,[0 L],'r'),...
+     title(['poutre appuyee - libre ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+   end
+ case 5, disp('poutre encastree - appuyee');
+  figure('Name','3 premiers modes de vibration d''une poutre encastree - appuyee',...
+  'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+  subplot(nmod+1,1,1),hold on, fplot('[tan(x) , tanh(x)]',[0,4*pi,-1.5,1.5],'b'),...
+  title('representation de l''equation implicite tan(x) = th(x)'), grid   
+  for imod = 1:nmod
+     a= imod*pi/2+.1;
+     b=(2*imod+1)*pi/2-.1;
+     lx = fzero(@(x) tan(x)-tanh(x),[a b])
+     anal=(lx^2)*sqrt(EI/roS/(L^4));
+     aa = (cos(lx)+cosh(lx))/(sin(lx)+sinh(lx));
+     mode = @(x) cos(lx*x/L)-cosh(lx*x/L)-aa*(sin(lx*x/L)-sinh(lx*x/L));
+     subplot(nmod+1,1,imod+1), fplot(mode,[0 L],'r'),...
+     title(['poutre encastree - appuyee ',num2str(imod),'ieme mode de pulsation ',num2str(anal),...
+     ' Hz ' ]), grid
+   end    
+end
+end
+disp('les autres cas seront simples a programmer en utilisant le cours');

Work/FV_4.m

@@ -0,0 +1,94 @@
+clc;clear all;close all;
+% comparaison des solutions numeriques et analytique de l'exercice FV4 
+% Barre bi-encastree
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+ES= 1 ; roS=1; L= 1;
+%----------------------------------------------------------
+%----- Methodes d'approximation sur la formulation forte
+%----- avec la fonction de forme x(L-x)
+omeg1=pi*sqrt(ES/roS/L/L);    % solution analytique
+taille = get(0,'ScreenSize'); 
+figure('Name','comparaison fonction de forme - premier mode',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+fplot(@(x)[sin(pi*x/L) , 4*x.*(L-x)],[0 L]);
+disp('---------------------------------------');
+disp('comparaison de la fonction de forme et du premier mode');
+disp('---------------------------------------');
+pause(2);
+close all;
+%Collocation
+f=@(x) x.*(L-x); % D2f=@(x) -2*x./x;
+disp('Methode de Collocation');
+xp=0.5;
+while xp > 0
+k = 2*ES ; m = roS*f(xp) ; sol = sqrt(k/m);
+disp('---------------------------------------');
+fprintf('Pour un point de collocation situe en %3.1f L \n',xp);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('pour quitter la methode de Collocation ne pas donner de nouvelle valeur');
+xp = input('donner une nouvelle position du point de collocation sur [0,1]: ');
+if isempty(xp) xp=0; end 
+end
+disp('---------------------------------------');
+disp('Methode de la valeur moyenne');
+k = 2*ES*L; m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode de galerkin');
+f2 = @(x) (x.*(L-x)).^2;
+k = 2*ES*quad(f,0,L); m = roS*quad(f2,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+close all;
+%----------------------------------------------------------
+%------ Methode des elements finis
+disp('---------------------------------------');
+disp('Methode des elements finis');
+ne=2;
+while ne > 1
+    K=zeros(ne-1);	M=zeros(ne-1);	 
+figure('Name',['comparaison EF - analytique pour un mailage de '...
+    ,num2str(ne),' elements'],...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+Le=L/ne;
+if ne <= 2 K(1,1)=2*ES/Le; M(1,1)=4*roS*Le/6;end
+for i=1:ne-2         
+    K(i,i)=2*ES/Le;  M(i,i)=4*roS*Le/6;
+    K(i,i+1)=-ES/Le; M(i,i+1)=roS*Le/6;
+    K(i+1,i)=-ES/Le; M(i+1,i)=roS*Le/6;
+    K(i+1,i+1)=2*ES/Le;  M(i+1,i+1)=4*roS*Le/6;
+end
+if ne >= 4 nmod=3; else nmod =ne-1; end
+[modes,omega] = eigs(K,M,ne-1,'sm'); 
+omeg = diag(sqrt(omega));
+disp('---------------------------------------');
+fprintf('Pour un maillage de %2d elements finis \n',ne);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',omeg(1),omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(omeg(1)-omeg1)/omeg1);
+
+for imod = 1:nmod
+    U= modes(:,imod)  ; maxi = max(abs(U));
+    if U(1)<0  U = -U ;end 
+    U = [ 0 ;U ;0]/maxi;
+% trace des modes 
+    for i=1:ne
+        x1=(i-1)*Le; x2=i*Le;
+        x=x1:Le/10:x2; y=U(i)+(U(i+1)-U(i))*(x-x1)/Le; 
+        subplot(nmod,1,imod),hold on, plot(x,y,'b');
+    end
+sol = omeg(imod); anal=imod*pi*sqrt(ES/roS/L/L);
+err = 100*abs(sol-anal)/anal;
+subplot(nmod,1,imod), fplot(@(x) sin(imod*pi*x/L),[0 L],'r'),...
+title([num2str(imod),'ieme mode : pulsation approchee ',num2str(sol),...
+  ' / analytique ',num2str(anal),' soit une erreur de ',num2str(err),...
+  ' % ' ]), legend('sol. analytique en rouge','Sol. EF en bleu'), grid
+end
+disp('pour quitter l''application ne pas donner de nouvelle valeur');
+ne = input('donner le nombre d''elements ne ? [>2]: ');
+if isempty(ne) ne=0; end 
+end
+close all;

Work/FV_5.m

@@ -0,0 +1,142 @@
+clc;clear all;close all;
+% comparaison des solutions numeriques et analytique de l'exercice FV5 
+% Mur climatise
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+T0= 50 ; phie= 10; ro=-10*pi; lamda=1; e=1;
+%----- solution analytique
+%ue = dsolve('D2u = -ro*sin(pi*t)/lamda','u(0)=T0','Du(e)=-phie/lamda');
+%ue = simple(ue)
+fa=@(x) T0+x*(-phie+ro/pi)+ro*sin(pi*x)/pi/pi ;
+%----------------------------------------------------------
+%----- Methode de Galerkin
+disp('---------------------------------------');
+disp('Methode de Galerkin');
+n = input('donner le degre n de l''approximation ? [2]: ');
+if isempty(n) n=2; end 
+disp('---------------------------------------');
+taille = get(0,'ScreenSize'); 
+figure('Name','comparaison Galerkin - analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+hold on,
+fplot(fa,[0 e],'r');
+K=[];F=[];
+for i=1:n
+    for j=1:n
+    F1 = @(x) i*j*x.^(i+j-2);
+    K(i,j)= quad(F1,0,e);
+    end
+    F2 = @(x) sin(x*pi/e).*x.^i;
+    F(i)= ro*quad(F2,0,e)-phie*e^i;
+end
+disp('---------------------------------------');
+fprintf('Methode de galerkin pour un polynome de degre %3i\n',n);
+disp('systeme matriciel');
+disp(K); disp(F);
+U = K\F' ;
+y =T0;
+x=0:e/100:e;
+for i=1:n   y=y+U(i)*x.^i;  end
+plot(x,y,'b');
+legend('sol. analytique en rouge','approx par Galerkin en bleu',2);
+grid
+%----------------------------------------------------------
+%------ Methode des elements finis
+disp('---------------------------------------');
+disp('Methode des elements finis');
+ne = input('donner le nombre d''elements ne ? [2]: ');
+if isempty(ne) ne=2; end 
+figure('Name','comparaison EF - analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+hold on, fplot(fa,[0 e],'r');
+Le=e/ne;
+K=[];F=[];
+for i=1:ne-1         
+    K(i,i)=2*lamda/Le;
+    K(i,i+1)=-lamda/Le;
+    K(i+1,i)=-lamda/Le;
+    F1 = @(x) (sin(pi*(i/ne+x/e)) - 2*sin(pi/(2*ne))*...
+    cos(pi*((i-0.5)/ne+x/e)).*x./Le);
+    F(i)= ro*quad(F1,0,Le);
+end
+F(1)=F(1)+T0*lamda/Le;
+K(ne,ne)=lamda/Le;
+F2 = @(x) sin(pi*((ne-1)/ne+x/e)).*x./Le;
+F(ne)= ro*quad(F2,0,Le)-phie;
+disp('---------------------------------------');
+fprintf('systeme matriciel pour %2d elements finis\n',ne);
+    disp(K); disp(F);
+U = K\F' ;
+x1=0;x2=Le;
+x=x1:Le/10:x2; y=T0+(U(1)-T0)*(x-x1)/Le; plot(x,y,'b');
+for i=2:ne
+    x1=(i-1)*Le; x2=i*Le;
+    x=x1:Le/10:x2; y=U(i-1)+(U(i)-U(i-1))*(x-x1)/Le; plot(x,y,'b');
+end
+legend('sol. analytique en rouge','Sol. EF en bleu',2);
+grid
+%----------------------------------------------------------
+%----- Methode de Collocation
+disp('---------------------------------------');
+disp('Methode de Collocation');
+n = input('donner le degre n de l''approximation ? [2]: ');
+if isempty(n) n=2; end 
+figure('Name','comparaison Collocation - analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+hold on,
+fplot(fa,[0 e],'r');
+K=[];F=[];
+F(1)=-phie ;
+for j=1:n  K(1,j)= j ; end
+for i=2:n
+    xi=e*(i-1)/n;
+    for j=1:n   K(i,j)= j*(j-1)*xi^(j-2); end
+    F(i)= -ro*sin(pi*xi/e)/lamda;
+end
+disp('---------------------------------------');
+fprintf('Methode de collocation pour un polynome de degre %3i\n',n);
+disp('systeme matriciel');
+disp(K); disp(F);
+U = K\F' ;
+y =T0;
+x=0:e/100:e;
+for i=1:n
+    y=y+U(i)*x.^i;
+end
+plot(x,y,'b');
+legend('sol. analytique en rouge','approx par Collocation en bleu',2);
+grid
+%----------------------------------------------------------
+%----- Methode de la valeur moyenne par sous domaine
+disp('---------------------------------------');
+disp('Methode de la valeur moyenne');
+n = input('donner le degre n de l''approximation ? [2]: ');
+if isempty(n) n=2; end 
+figure('Name','comparaison Valeur moyenne - analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+hold on,
+fplot(fa,[0 e],'r');
+nsd=n-1; Ld=e/nsd;
+K=[];F=[];
+F(1)=-phie ;
+for j=1:n  K(1,j)= j; end
+for i=2:n
+    e1=(i-2)*Ld; e2=e1+Ld;
+    for j=1:n   K(i,j)= j*(e2^(j-1)-e1^(j-1)); end
+    F2 = @(x) sin(x*pi/e);
+    F(i)= -ro*quad(F2,e1,e2)/lamda;
+end
+disp('---------------------------------------');
+fprintf('Methode de la valeur moyenne pour un polynome de degre %3i\n',n);
+disp('systeme matriciel');
+disp(K); disp(F);
+U = K\F' ;
+y =T0;
+x=0:e/100:e;
+for i=1:n
+    y=y+U(i)*x.^i;
+end
+plot(x,y,'b');
+legend('sol. analytique en rouge','approx par valeurs moyennes en bleu',2);
+grid

Work/FV_6.m

@@ -0,0 +1,85 @@
+clc;clear all;close all;
+% comparaison de la solution numerique et analytique de l'exercice FV6 
+% Probleme de convection-diffusion (stationaire)avec stabilisation SUPG
+% conditions de dirichlet nulles aux deux extremites
+%  H.Oudin 
+L=1;r=1;V=1; 
+Pe = input('donner le nombre de peclet du probleme ? [10]: ');
+if isempty(Pe) Pe=10; end   
+Dk=V*L/Pe;   % diffusivite du milieu en fonction de Pe
+  
+nelt = input('donner le nombre d''element en x (longueur) ? [3]: ');
+if isempty(nelt) nelt=3; end  
+n=nelt-1;
+Le=L/nelt;
+K=zeros(n); F=r*Le*ones(n,1);
+taille = get(0,'ScreenSize'); 
+figure('Name','comparaison sol EF - analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+hold on,
+%----- solution analytique il faut utiliser la solution ue dans fa1
+%ue = dsolve('-k*D2u + V*Du = r','u(0)=0','u(L)=0');
+%ue = simple(ue);
+fa1=@(t) r*(-exp(V/Dk*t)*L-t+t*exp(V/Dk*L)+L)/V/(-1+exp(V/Dk*L));
+fplot(fa1,[0 L],'r');
+%fa=@(x) x/L-(exp(Pe*x/L)-1)/(exp(Pe)-1);
+%fplot(fa,[0 L],'r'), 
+%----- solution numerique sans stabilisation
+Peh=V*Le/(2*Dk);
+for i=1:n-1          
+    K(i,i)=2/Peh;
+    K(i,i+1)=1-1/Peh;
+    K(i+1,i)=-1-1/Peh;
+end
+K(n,n)=2/Peh;    %disp(K); disp(F);
+K=V*K/2.;
+U = K\F ;
+x1=0;x2=Le;
+x=x1:Le/10:x2; y=U(1)*x/Le; plot(x,y,'b');
+for i=2:n
+    x1=(i-1)*Le; x2=i*Le;
+    x=x1:Le/10:x2; y=U(i-1)+(U(i)-U(i-1))*(x-x1)/Le; plot(x,y,'b');
+end
+x1=n*Le; x2=L;
+x=x1:Le/10:x2; y=U(n)-U(n)*(x-x1)/Le; plot(x,y,'b');
+legend('sol. analytique en rouge','Sol. EF en bleu',2);
+grid
+%----- solution numerique avec stabilisation SUPG
+
+
+beta = input('donner la valeur du coefficient de stabilisation choisie? [nominal]: ');
+if isempty(beta) beta=Le*(coth(Peh)-1/Peh)/(2*V)% valeur optimale de stabilisation 
+end   
+figure('Name','comparaison sol EF - analytique avec stabilisation',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+hold on,
+fa1=@(t) r*(-exp(V/Dk*t)*L-t+t*exp(V/Dk*L)+L)/V/(-1+exp(V/Dk*L));
+fplot(fa1,[0 L],'r');
+x1=0;x2=Le;
+x=x1:Le/10:x2; y=U(1)*x/Le; plot(x,y,'b');
+for i=2:n
+    x1=(i-1)*Le; x2=i*Le;
+    x=x1:Le/10:x2; y=U(i-1)+(U(i)-U(i-1))*(x-x1)/Le; plot(x,y,'b');
+end
+x1=n*Le; x2=L;
+x=x1:Le/10:x2; y=U(n)-U(n)*(x-x1)/Le; plot(x,y,'b');
+
+Ksup=zeros(n);
+for i=1:n-1          
+    Ksup(i,i)=2/Le;
+    Ksup(i,i+1)=-1/Le;
+    Ksup(i+1,i)=-1/Le;
+end
+Ksup(n,n)=2/Le;    %disp(K); disp(F);
+K = K + V*V*beta*Ksup;
+U1 = K\F ;
+x1=0;x2=Le;
+x=x1:Le/10:x2; y=U1(1)*x/Le; plot(x,y,'g');
+for i=2:n
+    x1=(i-1)*Le; x2=i*Le;
+    x=x1:Le/10:x2; y=U1(i-1)+(U1(i)-U1(i-1))*(x-x1)/Le; plot(x,y,'g');
+end
+x1=n*Le; x2=L;
+x=x1:Le/10:x2; y=U1(n)-U1(n)*(x-x1)/Le; plot(x,y,'g');
+legend('sol. analytique en rouge','Sol. EF en bleu','avec stabilisation en vert',2);
+grid

Work/NUM_4.m

@@ -0,0 +1,69 @@
+clc;clear all;close all;
+% Mapping d'un element Q4 etude de la transformation geometrique
+% Exercice NUM4 du chapitre NUM du site
+%  H.Oudin 
+taille = get(0,'ScreenSize');           
+%------ mappage de l'element de reference
+[S,T] = meshgrid(-1:0.2:1);
+[n,n] = size(S);
+%------ donnee de la geometrie de l'element reel (position des 4 noeuds)
+disp('---------------------------------------');
+disp('Mapping d''un element Q4 etude de la transformation geometrique');
+disp('cliquez les 4 noeuds de l''element reel dans cette figure');
+disp('---------------------------------------');
+figure('Name','donnees de l''element reel',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])    
+axis([0 1 0 1])
+grid on
+[Coord(:,1),Coord(:,2),button] = ginput(4); 
+close all;                                      
+%Coord=[0 0;8 0;8 2;4 7 ];% peut etre utilise comme donnee
+%------ mapping de l'element   
+for i=1:n
+    for j=1:n
+        s=S(j,i);t=T(j,i);
+        N = .25*[(1-s)*(1-t)  (1+s)*(1-t)  (1+s)*(1+t)  (1-s)*(1+t)];
+        %------ calcul du mapping de l'element reel 
+        Xr(j,i)=N*Coord(:,1);Yr(j,i)=N*Coord(:,2);
+     %------ calcul de detJ sur le mapping
+     dN = .25*[-(1-t)  (1-t) (1+t)  -(1+t)
+   		       -(1-s) -(1+s) (1+s)   (1-s)]; 
+     J = dN*Coord;
+     det(j,i) = J(1,1)*J(2,2)-J(1,2)*J(2,1);
+    end
+end
+C = max(det); maxi=max(C); D = min(det);mini=min(D);
+disp('Figures donnant l''evolution de la valeur de detJ');
+disp('Le mapping represente les lignes a s ou t constant');
+disp('Attention l''element reel est plan, ne vous laissez pas influencer par le trace');
+disp('-------------------------------------------------------------------------------');
+figure('Name','Valeur de detJ sur le mapping de l''element de reference',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])    
+subplot(1,2,1), surf(S,T,det),title('element de reference'),axis equal ...
+    , colormap jet, colorbar, view([0 0 1])
+subplot(1,2,2), surf(Xr,Yr,det),title('element reel'), axis equal, ...
+      view([0 0 1])
+% calcul de la surface de l'element par integration numerique
+npg = 4;          %----- integration a 4 points de Gauss
+wg = [1,1,1,1];   %----- poids et position
+c = 1/sqrt(3); posg = [ -c -c ; c -c ; c  c ; -c  c ];
+aire=0;
+for ipg=1:npg       %----- boucle d'integration
+   s = posg(ipg,1); t = posg(ipg,2); poids = wg(ipg);
+   dN = .25*[-(1-t)  (1-t) (1+t)  -(1+t)
+   		     -(1-s) -(1+s) (1+s)   (1-s)]; 
+   J = dN*Coord;
+   detj = J(1,1)*J(2,2)-J(1,2)*J(2,1);
+   aire=aire+detj*poids;  
+end
+Coord 
+if sign(maxi) == sign(mini) disp('---------------------------------');
+    disp('detJ ne s''annule pas sur l''element reel');
+    disp('la transformation geometrique est bijective ');
+    fprintf('Par integration numerique la surface de l''element reel est : %6.4f\n',aire);
+else disp('---------------------------------');
+    disp('detJ change de signe sur l''element');
+    disp('la transformation geometrique n''est pas bijective ');
+    fprintf('La valeur mini du determinant est : %6.4f  et le max : %6.4f \n',mini,maxi);
+    fprintf('L''integration numerique donne une surface calculee de %6.4f\n',aire);
+end

Work/barre_approx_cours.m

@@ -0,0 +1,94 @@
+clc;clear all;close all;
+% comparaison des solutions numeriques et analytique de l'exercice de cours 
+% Vibration d'une barre encastree - libre
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+ES= 1 ; roS=1; L= 1;
+%----------------------------------------------------------
+%----- Methodes d'approximation sur la formulation forte
+%----- avec la fonction de forme x(2x-L)
+disp('---------------------------------------');
+disp('Vibration d''une barre encastree - libre');
+disp('---------------------------------------');
+disp('Methodes d''approximation avec la fonction de forme x(2x-L)');
+disp('Utilise la formulation forte du probleme');
+omeg1=pi*sqrt(ES/roS/L/L)/2;    % solution analytique
+taille = get(0,'ScreenSize'); 
+figure('Name','comparaison de la fonction de forme et du premier mode',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+fplot(@(x)[sin(pi*x/2/L) , x.*(2*L-x)],[0 L]);
+%Collocation
+f=@(x) x.*(2*L-x); % D2f=@(x) -2*x./x;
+disp('---------------------------------------');
+disp('Methode de Collocation');
+xp=0.5;
+while xp > 0
+  k = 2*ES ; m = roS*f(xp) ; sol = sqrt(k/m);
+  disp('---------------------------------------');
+  fprintf('Pour un point de collocation situe en %3.1f L \n',xp);
+  fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+  fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+  disp('---------------------------------------');
+  disp('pour quitter la methode de Collocation ne pas donner de nouvelle valeur');
+  xp = input('donner une nouvelle position du point de collocation sur [0,1]: ');
+  if isempty(xp) xp=0; end 
+end
+disp('---------------------------------------');
+disp('Methode de la valeur moyenne');
+k = 2*ES*L; m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode de galerkin');
+f2 = @(x) (x.*(2*L-x)).^2;
+k = 2*ES*quad(f,0,L); m = roS*quad(f2,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+close all;
+%----------------------------------------------------------
+%----- Méthode de Galerkin sur la formulation faible pour un polynome de degre n
+disp('---------------------------------------------------------------------------');
+disp('---------------------------------------------------------------------------');
+n=2;
+fprintf('Methode de Galerkin sur la formulation faible pour un polynome de degre %3i\n',n);
+while n > 1
+K=zeros(n);	M=zeros(n);	
+for i=1:n
+    for j=1:n
+     F1 = @(x) i*j*x.^(i+j-2); K(i,j)= ES*quad(F1,0,L);
+     F2 = @(x) x.^(i+j); M(i,j)= roS*quad(F2,0,L);
+    end
+end
+%disp('systeme matriciel');K M
+if n >= 4 nmod=3; else nmod =n; end
+[modes,omega] = eigs(K,M,n,'sm'); 
+omeg = diag(sqrt(omega));
+%modes
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',omeg(1),omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(omeg(1)-omeg1)/omeg1);
+
+figure('Name','comparaison Galerkin - analytique',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+
+for imod = 1:nmod
+  U= modes(:,imod)  ; %
+  %if U(1)<0  U = -U ;end 
+  y =0;   x=0:L/100:L;
+  for i=1:n  y=y+U(i)*x.^i;  end
+  maxi = max(abs(y)); y=y/maxi;
+  subplot(nmod,1,imod),hold on, plot(x,y,'b'),plot(x,-y,'b')
+  sol = omeg(imod); anal=(2*imod-1)*pi*sqrt(ES/roS/L/L)/2;
+  err = 100*abs(sol-anal)/anal;
+  subplot(nmod,1,imod), fplot(@(x) -sin((2*imod-1)*pi*x/L/2),[0 L],'r'),
+  subplot(nmod,1,imod), fplot(@(x) sin((2*imod-1)*pi*x/L/2),[0 L],'r'),...
+  title([num2str(imod),'ieme mode : pulsation approchee ',num2str(sol),...
+  ' / analytique ',num2str(anal),' soit une erreur de ',num2str(err),...
+  ' % ' ]), legend('sol. analytique en rouge','Sol. approchee en bleu'), grid
+end
+disp('---------------------------------------');
+disp('pour quitter l''application ne pas donner de nouvelle valeur');
+n = input('donner le degre n de l''approximation ? [2]: ');
+if isempty(n) n=0; end 
+end
+close all;

Work/barre_approx_exo3.m

@@ -0,0 +1,74 @@
+clc;clear all;close all;
+%----- Methodes d'approximation sur la formulation forte
+% comparaison des solutions numeriques et analytique de l'exercice de cours 03
+% Vibration d'une barre encastree - ressort
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+ES= 1 ; roS=1; L= 1;
+alpa = input('donner la valeur du rapport ES/kl ? [10]: ');
+if isempty(alpa) alpa=10; end
+a= pi/2+0.01;                    % solution analytique
+b=3*pi/2-0.01;
+lx = fzero(@(x) tan(x)+alpa*x,[a b]);
+omeg1=lx*sqrt(ES/roS)/L
+%----------------------------------------------------------
+disp('Methodes d''approximation sur la formulation forte');
+disp('avec une fonction de comparaison polynomiale de degre 2');
+disp('---------------------------------------');
+beta = (alpa+1)/(2*alpa+1)
+f=@(x) x.*(1-beta*x/L);
+disp('Methode de Collocation');
+xp=0.5;
+while xp > 0
+  k = 2*beta*ES/L ; m = roS*f(xp*L) ; sol = sqrt(k/m);
+  disp('---------------------------------------');
+  fprintf('Pour un point de collocation situe en %3.1f L \n',xp);
+  fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+  fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+  disp('pour quitter la methode de Collocation ne pas donner de nouvelle valeur');
+  xp = input('donner une nouvelle position du point de collocation sur [0,1]: ');
+  if isempty(xp) xp=0; end 
+end
+disp('---------------------------------------');
+disp('Methode de la valeur moyenne');
+k = 2*beta*ES ; m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode de galerkin');
+f2 = @(x) (x.*(1-beta*x/L)).^2;
+k = 2*beta*ES*quad(f,0,L)/L ; m = roS*quad(f2,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+%----------------------------------------------------------
+disp('------------------------------------------');
+disp('amelioration de la fonction de comparaison de forme sinusoidale ');
+disp('------------------------------------------');
+beta=-(1+alpa*pi/2)
+f=@(x) 1 - cos(pi*x/2/L)+beta*sin(pi*x/2/L);
+D2f=@(x) (pi/2/L)^2*(cos(pi*x/2/L)-beta*sin(pi*x/2/L));
+disp('Methode de Collocation');
+xp=0.5;
+while xp > 0
+  k = -ES*D2f(xp*L) ; m = roS*f(xp*L) ; sol = sqrt(k/m);
+  disp('---------------------------------------');
+  fprintf('Pour un point de collocation situe en %3.1f L \n',xp);
+  fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+  fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+  disp('pour quitter la methode de Collocation ne pas donner de nouvelle valeur');
+  xp = input('donner une nouvelle position du point de collocation sur [0,1]: ');
+  if isempty(xp) xp=0; end 
+end
+disp('---------------------------------------');
+disp('Methode de la valeur moyenne');
+k = -ES*quad(D2f,0,L) ; m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode de galerkin');
+f1 = @(x) (1 - cos(pi*x/2/L)+beta*sin(pi*x/2/L)).*(pi/2/L)^2*(cos(pi*x/2/L)-beta*sin(pi*x/2/L));
+f2 = @(x) (1 - cos(pi*x/2/L)+beta*sin(pi*x/2/L)).*(1 - cos(pi*x/2/L)+beta*sin(pi*x/2/L));
+k = -ES*quad(f1,0,L) ; m = roS*quad(f2,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);

Work/barre_approx_exo4.m

@@ -0,0 +1,47 @@
+clc;clear all;close all;
+%----- Methodes d'approximation sur la formulation faible
+% comparaison des solutions numeriques et analytique de l'exercice de cours 04
+% Vibration d'une barre encastree - ressort
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+ES= 1 ; roS=1; L= 1;
+alpa = input('donner la valeur du rapport ES/kl ? [10]: ');
+if isempty(alpa) alpa=10; end 
+%taille = get(0,'ScreenSize'); 
+%figure('Name','solution analytique : 3 premiers modes de vibration de la barre',...
+%        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])                
+for imod = 1:3    % solution analytique
+  a= (2*imod-1)*pi/2+0.01;
+  b=(2*imod+1)*pi/2-0.01;
+  lx = fzero(@(x) tan(x)+alpa*x,[a b]);
+  anal(imod)=lx*sqrt(ES/roS/L/L);
+  %subplot(3,1,imod), fplot(@(x) sin(lx*x/L),[0 L],'r'),...
+  %title(['barre encastree - ressort : ',num2str(imod),...
+  %' ieme mode de pulsation ',num2str(anal(imod)),' Hz ' ]), grid
+end
+n=2;
+while n > 1
+K=zeros(n);	M=zeros(n);	            
+for i=1:n 
+  for j=1:n  
+    F1 = @(x) i*j*x.^(i+j-2);F2 = @(x) x.^(i+j);
+    K(i,j)= ES*(quad(F1,0,L)+(L^(i+j-1))/alpa); 
+    M(i,j)= roS*quad(F2,0,L); 
+  end    
+end
+if n >= 4 nmod=3; else nmod =n; end
+[modes,omega] = eigs(K,M,nmod,'sm'); 
+omeg = diag(sqrt(omega));
+disp('--------------------------------------------------');
+fprintf('Methode de galerkin pour un polynome de degre %2d \n',n);
+for imod = 1:nmod
+  fprintf('%2d pulsation propre : Approximation : %6.3f / analytique %6.3f\n',imod,omeg(imod),anal(imod));
+  fprintf('Soit une erreur de %2.1f %%\n',100*abs(omeg(1)-anal(1))/anal(1));
+end
+disp('---------------------------------------');
+disp('pour quitter l''application ne pas donner de nouvelle valeur');
+n = input('donner le degre du polynome ? [>2]: ');
+if isempty(n) n=0; end 
+end
+close all;

Work/barre_vibra2.m

@@ -0,0 +1,57 @@
+clc;clear all;close all;
+% Methode d'approximation appliquee aux vibrations des barres
+% comparaison des solutions numeriques et analytique de l'exercice vibra2 
+% Barre bi-encastree
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+ES= 1 ; roS=1; L= 1;
+%----------------------------------------------------------
+%----- Methodes d'approximation sur la formulation forte
+%----- avec la fonction de forme x(L-x)
+omeg1=pi*sqrt(ES/roS/L/L);    % solution analytique
+taille = get(0,'ScreenSize'); 
+%figure('Name','comparaison fonction de forme - premier mode',...
+%        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2])            
+%fplot(@(x)[sin(pi*x/L) , 4*x.*(L-x)],[0 L]);
+disp('Methodes d''approximation sur la formulation forte pour une barre bi-encastree');
+disp('------------------------------------------------------------------------------');
+disp('Methode de la valeur moyenne');
+f=@(x) x.*(L-x);
+k = 2*ES*L; m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode des residus ponderes avec p(x)=x');
+f=@(x) x.^2*(L-x);
+k = 2*ES*quad(@(x) x,0,L); m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode des residus ponderes avec p(x)=x^2');
+f=@(x) x.^3*(L-x);
+k = 2*ES*quad(@(x) x^2,0,L); m = roS*quad(f,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+disp('---------------------------------------');
+disp('Methode de Collocation');
+f=@(x) x.*(L-x); xp=0.5;
+while xp > 0
+  k = 2*ES ; m = roS*f(xp) ; sol = sqrt(k/m);
+  fprintf('Pour un point de collocation situe en %3.1f L \n',xp);
+  fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+  fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+  disp('---------------------------------------');
+  disp('pour quitter la methode de Collocation ne pas donner de nouvelle valeur');
+  xp = input('donner une nouvelle position du point de collocation sur [0,1]: ');
+  if isempty(xp) xp=0; end 
+end
+disp('---------------------------------------');
+disp('Methode de galerkin');
+f=@(x) x.*(L-x); f2 = @(x) (x.*(L-x)).^2;
+k = 2*ES*quad(f,0,L); m = roS*quad(f2,0,L); sol = sqrt(k/m);
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+close all;
+%----------------------------------------------------------
+close all;

Work/poutre_approx_exo9.m

@@ -0,0 +1,76 @@
+clc;clear all;close all;
+% comparaison des solutions numeriques et analytique de l'exercice de cours exo9
+% Vibration d'une poutre encastree - appuyee
+%  H.Oudin 
+%----------------------------------------------------------
+% donnees du probleme
+EI= 1 ; roS=1; L= 1;
+
+% solution analytique pour le premier mode
+lx1 = fzero(@(x) tan(x)-tanh(x),[pi/2+.1 3*pi/2-.1]);
+omeg1=(lx1^2)*sqrt(EI/roS/(L^4));
+aa = (cos(lx1)+cosh(lx1))/(sin(lx1)+sinh(lx1));
+mode1 = @(x) -cos(lx1*x/L)+cosh(lx1*x/L)+aa*(sin(lx1*x/L)-sinh(lx1*x/L));
+disp('---------------------------------------');
+disp('Methode de Galerkin a un parametre');
+%  Methode de Galerkin a un parametre (formulation forte)
+%  avec les fonctions de forme w1 = x^2(L-x)(3L-2x) et w2 = x^3(L-x)(4L-3x)
+w1 = @(x) ((x.^2).*(L-x)).*(3*L-2*x); 
+k(1,1) = 48*EI*quad(w1,0,L);
+w1w1 = @(x) (((x.^2).*(L-x)).*(3*L-2*x)).*(((x.^2).*(L-x)).*(3*L-2*x));
+m11 = roS*quad(w1w1,0,L);
+w1 = @(x) (((x.^2).*(L-x)).*(3*L-2*x))/sqrt(m11);
+m(1,1)= m11;
+sol = sqrt(k(1,1)/m11);
+disp('---------------------------------------');
+fprintf('Pour W1 l''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+%sol*sol
+w2 = @(x) ((x.^3).*(L-x)).*(4*L-3*x);
+k(2,2) = EI*quad(@(x) (((x.^3).*(L-x)).*(4*L-3*x)).*(360*x-168*L),0,L);
+w2w2 = @(x) (((x.^3).*(L-x)).*(4*L-3*x)).*(((x.^3).*(L-x)).*(4*L-3*x));
+m22 = roS*quad(w2w2,0,L); m(2,2) = m22;
+w2 = @(x) (((x.^3).*(L-x)).*(4*L-3*x))/sqrt(m22);
+sol = sqrt(k(2,2)/m22);
+disp('---------------------------------------');
+fprintf('Pour W2 l''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',sol,omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(sol-omeg1)/omeg1);
+%sol*sol
+taille = get(0,'ScreenSize'); 
+figure('Name','comparaison des fonctions de forme et du premier mode',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+hold on
+fplot(mode1,[0 L],'r'); fplot(w1,[0 L],'b');fplot(w2,[0 L],'g'),...
+,legend('sol. analytique en rouge , fonction W1 en bleu ,et W2 en vert '), grid
+pause(1);
+disp('---------------------------------------');
+disp('---------------------------------------');
+disp('Methode de Galerkin a deux parametres');
+disp('---------------------------------------');
+% La solution analytique pour le second mode
+lx2 = fzero(@(x) tan(x)-tanh(x),[3*pi/2+.1 5*pi/2-.1]);
+omeg2=(lx2^2)*sqrt(EI/roS/(L^4));
+aa = (cos(lx2)+cosh(lx2))/(sin(lx2)+sinh(lx2));
+mode2 = @(x) -cos(lx2*x/L)+cosh(lx2*x/L)+aa*(sin(lx2*x/L)-sinh(lx2*x/L));
+%----- Méthode de Galerkin a 2 parametres
+k(1,2) = EI*quad(@(x) (((x.^2).*(L-x)).*(3*L-2*x)).*(360*x-168*L),0,L);
+k(2,1) = 48*EI*quad(@(x) ((x.^3).*(L-x)).*(4*L-3*x),0,L);
+%disp('matrice K');k
+w1w2 = @(x) (((x.^2).*(L-x)).*(3*L-2*x)).*(((x.^3).*(L-x)).*(4*L-3*x));
+m(1,2) = roS*quad(w1w2,0,L); m(2,1)= m(1,2);
+%disp('matrice M'); m
+[modes,omega] = eigs(k,m,2,'sm'); 
+omeg = diag(sqrt(omega));
+fprintf('L''approximation de la premiere pulsation propre est %6.3f / analytique %6.3f\n',omeg(1),omeg1);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(omeg(1)-omeg1)/omeg1);
+fprintf('L''approximation de la seconde pulsation propre est %6.3f / analytique %6.3f\n',omeg(2),omeg2);
+fprintf('Soit une erreur de %2.1f %%\n',100*abs(omeg(2)-omeg2)/omeg2);
+Z1 = @(x) -5.58*modes(1,1)*(((x.^2).*(L-x)).*(3*L-2*x))-...
+5.58*modes(2,1)*((x.^3).*(L-x)).*(4*L-3*x);
+Z2 = @(x) -33.5*modes(1,2)*(((x.^2).*(L-x)).*(3*L-2*x))-...
+33.5*modes(2,2)*((x.^3).*(L-x)).*(4*L-3*x);
+figure('Name','comparaison Analytique / Galerkin pour les deux modes ',...
+        'Position',[taille(3)/2.01 taille(4)/2.6 taille(3)/2 taille(4)/2]) 
+hold on
+fplot(mode1,[0 L],'r');fplot(mode2,[0 L],'r'); fplot(Z1,[0 L],'b');fplot(Z2,[0 L],'b'),...
+,legend('sol. analytique en rouge , Approximation en bleu '), grid