Logiciel d'initiation 
à la méthode des éléments finis

(Gilbert Touzot, février 1997)
_________________________

Introduction

Ce document met en œuvre la méthode des éléments finis, dune manière cohérente avec le livre "Une présentation de la méthode des éléments finis" de G.Dhatt et G.Touzot, édité par Maloine.
La méthode des éléments finis est une méthode de discrétisation des équations aux dérivées partielle, très utilisée dans la simulation en sciences de l'ingénieur. Elle s'applique en particulier à la mécanique des solides et des fluides, à la thermique, à l'électromagnétique et à l'électrodynamique, à l'hydraulique,...

Très brièvement la méthode des éléments consiste en les étapes suivantes :

Etape 1 : Entrée des données qui définissent le problème traité

- Définition du domaine de définition (1, 2 ou 3 dimensions) des équations par la spécification :
         - de points de discrétisation (ou nœuds) en lesquels seront définies les grandeurs inconnues
         - de sous domaines de formes simples (ou éléments) telles que des segments de droite à 1 dimension, des 
              triangles à 2 dimensions, des tétraèdres ou hexaèdres à 3 dimensions. 
         Chaque élément est défini par une liste de nœuds.
   - Définition des propriétés physiques des matériaux qui consituent les milieux étudiés
   - Définition des sollicitations externes auxquels sont soumis les milieux étudiés : forces, sources de chaleur, ...
   - Définition des "conditions aux limites" imposées sur la portion du contour du domaine étudié  sur laquelle la
          valeur des inconnues est imposée.

Etape 2 : Construction des matrices et vecteurs qui définissent le problème 
               discrétisé

Construction, selon le type de problème posé (stationnaire ou non, premier ou second ordre 
    par rapport au temps) de :
    
               - matrice de rigidité
               - matrice de masse
               - matrice d'amortissement
               - vecteur second membre
               - ...
   Cette étape fait intervenir les équations aux dérivées partielles qui régissent le système  physique
   étudié. Ces équations sont exprimées sous une forme faible (forme variationnelle) et discrétisée, 
   grâce à la technique"d'interpolation par éléments".

Etape 3 : Résolution du problème discrétisé par rapport aux variables d'espace

Selon le problème considéré, cette étape peut consister en :
    
             - la résolution d'un système d'équations algébriques linéaires, par une méthode directe
               ou itérative
             - calcul des valeurs et vecteurs propres d'un système
             - la résolution itérative d'un système d'équations algébriques non linéaires par une 
               méthode du type Newton-Raphson
             - l'intégration par rapport au temps d'un système d'équations différentielles par rapport 
               au temps, du premier ou du second ordre

Etape 4 : Post-traitement

Cette étape peut consister en :
    
             - le calcul, à partir de la solution obtenue, de variables secondaires, en général liées
               aux dérivées de la solution par rapport aux variables d'espace : déformations, 
               contraintes, gradients de pression ou de températures, ...
             - le calcul des résidus de déséquilibre nodaux, qui permettent d'évaluer la qualité de la 
               solution, et d'obtenir les réactions correspondant à des degrés de liberté imposés par des
               conditions aux limites
             - la présentation graphique des différents résultats.

Exemples :

- élément de type 1 : barre élastique à 1 dimension et 2 nœuds
             - élément de type 2 : barre élastique à 3 dimension et 2 nœuds
             - élément de type 3 : triangle élastique à 2 dimensions et 3 nœuds

Remarques :

- pour exécuter l'exemple numéro "n", il suffit d'imposer "nexemple = n"
             - les informations  qui constituent des données d'entrée sont commentées en rouge

Etape 1 : entrée des données qui définissent le problème étudié

Entrée des informations sur les noeuds

(* éviter les messages d'erreur inutiles  *) 
Off[General::spell1];

(* entrer le numéro de l'exemple à exécuter (entre 1 et 3) 
             1 : élément de type 1 : barre élastique à 1 dimension et 2 nœuds
             2 : élément de type 2 : barre élastique à 3 dimension et 2 nœuds
             3 : élément de type 3 : triangle Laplace à 2 dimensions et 3 nœuds *)
nexemple = 1;

(* entrer l'indice d'impression (entre 0 et 3)
           0 : pas d'impression
           1 : solution
           2 : vecteurs
           3 : matrices et vecteurs  *)
print = 3;

(* entrer les coordonnées des noeuds *)
Switch[nexemple,
                      (* modèle à 1 dimension *)
    1,  xn = {{0.},{0.5},{1.},{1.5}},
        (* modèle à 3 dimensions *)
    2, xn = {{0.,0.,0.},{1.,0.,0.},{0.,1.,0.},{1.,1.,0.},
             {0.5,0.5,1.}},
        (* modèle à 2 dimensions *)
    3, xn = {{0.,0.},{1.,0.},{0.,1.},{1.,1.},{0.5,0.5}} ];

(* entrer le nombre (constant) d'inconnues par noeud *)
Switch[nexemple,
    1, ndln = 1,
    2, ndln = 3,
    3, ndln = 1  ];



(* impression des coordonnées des nœuds *)
If[print > 1,
Print["coordonnées des nœuds :"];
Print[xn],];

(* calculs auxiliaires *)
nnt  = Dimensions[xn][[1]];(* nombre de noeuds *)
ndim = Dimensions[xn][[2]];(* nombre de dimensions d'espace *)
neq  = ndln * nnt;         (* nombre total d'inconnues *)

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Entrée des informations sur les éléments

(* entrer le type d'élément utilisé (supposé ici constant dans tout le modèle) *)
Switch[nexemple,
    1, ntpe = 1,
    2, ntpe = 2,
    3, ntpe = 3  ];

(* entrer les numéros des noeuds de chaque élément (connectivité )*)
Switch[nexemple,
           1, conn = {{1,2},{2,3},{3,4}},
    2, conn = {{1,5},{2,5},{3,5},{4,5}},
    3, conn = {{1,2,5},{1,5,3},{3,5,4},{2,4,5}}];

(* impression de la table de connectivité des éléments *)
If[print > 1,
Print["table de connectivité :"];
Print[conn],];

(* calculs auxiliaires *)
nel  = Dimensions[conn][[1]]; (* nombre d'éléments *)
nnel = Dimensions[conn][[2]]; (* nombre de noeuds (constant) de chaque élément *)

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Entrée des informations sur les propriétés physiques des matériaux

entrer les propriétés élémentaires, qui dépendent du type d'élément utilisé :

       - élément 1 (barre ou Laplace à 2 nœuds et 1 dimension) :
                                             module d'young (ou conductivité thermique)
                                             section de la barre
                                             masse volumique du  matériau
                                             force (ou source) par unité de longueur, suivant l'axe x local
                                                                                
     - élément 2 (barre ou Laplace à 2 nœuds et 3 dimensions) :
                                            module d'young (ou conductivité thermique)
                                            section de la barre
                                            masse volumique du  matériau
                                            force (ou source) par unité de longueur, suivant l'axe x local
     - élément 3 (triangle Laplace à 2 dimensions) : 
                                            conductivité thermique
                                            épaisseur
                                            masse volumique du  matériau
                                            source par unité de surface

(* entrer la table des propriétés élémentaires, ici identique pour tous les éléments *)
 Switch[nexemple,
           1, prel = {{1., 1., 1., 0.}},
    2, prel = {{1., 1., 1., 0.}},
    3, prel = {{1., 1., 1., 1.}}];

(* impressionde la table des propriétés élémentaires *)
If[print > 1,
Print[" table des propriétés élémentaires :"];
Print[prel],];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Entrée des informations sur les sollicitations concentrées

(* entrer pour, chaque sollicitation concentrée un triplet :
           {numéro du noeud sollicité, 
             numéro du degré de liberté (ddl) sollicité,
             valeur de la sollicitation concentrée appliquée au degré deliberté}    *)
  Switch[nexemple,
            1, sollicit = {{4,1,1.}},
     2, sollicit = {{5,2,-1.}},
     3, sollicit = {{2,1,1.}, {4,1,1.}}];

(* impression de la table des sollicitations *)
If[print > 1,
Print["table des sollicitations :"];
Print[sollicit],];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Entrée des informations sur les conditions aux limites

(* entrer, chaque condition à la limite, un triplet :
        {numéro du noeud à degré de liberté (ddl) imposé, 
          numéro du degré de liberté imposé,
          valeur imposée au degré de liberté}    *)
 Switch[nexemple,
           1,  condlim = {{1,1,1.}},
    2,  condlim = {{1,1,0.},{1,2,0.},{1,3,0.},
           {2,1,0.},{2,2,0.},{2,2,0.},
           {3,1,0.},{3,2,0.},{3,3,0.},
           {4,1,0.},{4,2,0.},{4,3,0.}},
    3,  condlim = {{1,1,0.}, {3,1,0.}}];

(* impression de la table des conditions aux limites *)
If[print > 1,
Print["table des conditions aux limites :"];
Print[condlim],];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Etape 2 : construction des matrices et vecteurs globaux

Fonctions de calcul spécifiques de divers types d'éléments :
          - matrices de rigidité élémentaires
          - post-calcul des gradients
          - tracé des résultats

Elément de type 1 : Laplace (ou barre) linéaire à 1 dimensions et à 2 noeuds

(* matrice de rigidité élémentaire *)
Clear[rgdt1];
rgdt1[prelel_,xnel_]:=Module[{},

    long = Abs[xnel[[2,1 ]] - xnel[[1,1 ]]];
    coef = prelel[[1]] * prelel[[2]] / long;
    ke = coef * {{1.,-1.},{-1.,1.}} ;
    Return[ke] ];
    
(* matrice de rigidité tangente élémentaire *)
Clear[rtgt1];
rtgt1[prelel_,xnel_]:=Module[{},
   (* non implantée *)
    Return[ke] ];
    
(* matrice de masse élémentaire *)
Clear[mass1];
mass1[prelel_,xnel_]:=Module[{},

    long = Abs[xnel[[2,1 ]] - xnel[[1,1 ]]];
    coef = prelel[[2]] * prelel[[3]] * long / 6.;
    me = coef * {{2.,1.},{1.,2.}} ;
    Return[me] ];
    
(* forces nodales dues à des sollicitations réparties *)
Clear[solr1];
solr1[prelel_,xnel_]:=Module[{},

    long = Abs[xnel[[2,1 ]] - xnel[[1,1 ]]];
    coef = prelel[[4]] * long / 2.;
    fe = coef * {1.,1.} ;
    Return[fe] ];
    
(*forces nodales dues à la déformation, qui fournissent le résidu  d'équilibre par assemblage*)
Clear[res1];
res1[prelel_,xnel_,unel_]:=Module[{},

    (* longueur de l'élément *)
    long = Abs[xnel[[2,1]] - xnel[[1,1]]];
    (* gradient, ou déformation axiale *)
    eps = (unel[[2]] - unel[[1]]) / long; 
    coef = prelel[[1]] * prelel[[2]];
    re = coef * eps * {-1,1};
    (* prise en compte des sollicitations réparties *)
           If[prelel[[4]] != 0.0, 
         cc = prelel[[4]] * long / 2.;
         re = re - cc * {1,1}, ];
    (* énergie de déformation *)
    energ = coef * eps * eps * long / 2.;
    (* impression des contraintes et déformations de l'élément *)
    Return[{energ,re}] ];
    
(* calcul et impression des variables de l'élément liées aux gradients : 
                   déformation, effort et contrainte axiaux *)
Clear[grad1];
grad1[prelel_,xnel_,unel_]:=Module[{},

    (* longueur de l'élément *)
    long = Abs[xnel[[2,1]] - xnel[[1,1]]];
    (* gradient, ou déformation axiale *)
    eps = (unel[[2]] - unel[[1]]) / long;
    coef = prelel[[1]] * prelel[[2]];
    (* forces axiale *)
    fax = coef * eps;
    (* contrainte *)
    sigma = fax / prelel[[2]];
    (* impression *)
    Print[" élément : ",iel];
    Print["      déformation  = ",eps];
    Print["      effort axial = ",fax];
    Print["      contrainte   = ",sigma];
    Return[{eps,fax,sigma}]  ];

(* tracé de l'élément avec déformée éventuelle *)
Clear[tracem1];
tracem1[xnel_,unel_,plot_]:=Module[{},

    (* coordonnées des nœuds de l'élément *)
    x1 = xnel[[1,1]]; x2 = xnel[[2,1]];
    (* tracé en "nplot"sous-éléments *)
    plot1=plot;
    plot1=Append[plot1,{Hue[0.3,1,1],Thickness[0.01],
          Line[{{x1,0.},{x2,0.}}]} ];  
    Return[plot1] ];
    
(* tracé d'une grandeur définie aux noeuds par interpolation linéaire *)
Clear[tracei1];
tracei1[xnel_,unel_,ampli_,umin_,umax_,ncoupe_,
        uplot_,plot_]:=Module[{},

    (* coordonnées des nœuds de l'élément *)
    x1 = xnel[[1,1]] + ampli * unel[[1,1]];
    x2 = xnel[[2,1]] + ampli * unel[[2,1]];
    dx = (x2 - x1)/ncoupe;
    (* valeurs aux nœuds de l'élément de la fonction à tracer  *)
    upl1 = uplot[[1]]; upl2 = uplot[[2]];
    dupl = (upl2 - upl1)/ncoupe;
    
    (* tracé en "nplot"sous-éléments *)
    plot1=plot;
    x = x1; upl = upl1 + dupl/2.;
    hmax = 0.9; (* coefficient de réglage de la couleur, entre 0 et 1 *)
    Do[ 
      If[ umax!=umin,
           h = 1. - hmax + hmax *(upl - umin) / (umax - umin),
           h = 1.];
      plot1=Append[plot1,{Hue[h,1,1],Thickness[0.01],
                              Line[{{x,0.},{x + dx,0.}}]} ];  
     x = x +dx; upl = upl + dupl,
      {i,1,ncoupe}];
    (* tracé *) 
     Return[plot1] ];

Elément de type 2 : Laplace (ou barre) linéaire à 3 dimensions et à 2 noeuds

(* matrice de rigidité élémentaire *)
Clear[rgdt2];
rgdt2[prelel_,xnel_]:=Module[{},

    dx = xnel[[2,1]] - xnel[[1,1]];
    dy = xnel[[2,2]] - xnel[[1,2]];
    dz = xnel[[2,3]] - xnel[[1,3]];
    long = Sqrt[dx^2 + dy^2 + dz^2];
    transf = {dx, dy, dz} / long;
    coef = prelel[[1]] * prelel[[2]] / long;
    c = Outer[Times, transf, transf] * coef;
    ke = Table[0.,{i,1,6},{j,1,6}];
    Do[
      Do[
          ke[[i    ,j    ]] =   c[[i,j]];
          ke[[i    ,j + 3]] = - c[[i,j]];
          ke[[i + 3,    j]] = - c[[i,j]];
          ke[[i + 3,j + 3]] =   c[[i,j]],
        {j,1,3}],
      {i,1,3}];
    Return[ke] ];
    
(* matrice de rigidité tangente élémentaire *)
Clear[rtgt2];
rtgt2[prelel_,xnel_]:=Module[{},
   (* non implantée *)
    Return[ke] ];
    
(* matrice de masse élémentaire *)
Clear[mass2];
mass2[prelel_,xnel_]:=Module[{},

    dx = xnel[[2,1]] - xnel[[1,1]];
    dy = xnel[[2,2]] - xnel[[1,2]];
    dz = xnel[[2,3]] - xnel[[1,3]];
    long = Sqrt[dx^2 + dy^2 + dz^2];
    coef = prelel[[2]] * prelel[[3]] * long /6.;
    me = Table[0.,{i,1,6},{j,1,6}];
    Do[
      Do[
          me[[i    ,j    ]] = 2. * coef;
          me[[i    ,j + 3]] = coef;
          me[[i + 3,    j]] = coef;
          me[[i + 3,j + 3]] = 2. * coef,
        {j,1,3}],
      {i,1,3}];
    Return[me] ];
    
(* forces nodales dues à des sollicitations réparties *)
Clear[solr2];
solr2[prelel_,xnel_]:=Module[{},

    dx = xnel[[2,1]] - xnel[[1,1]];
    dy = xnel[[2,2]] - xnel[[1,2]];
    dz = xnel[[2,3]] - xnel[[1,3]];
    long = Sqrt[dx^2 + dy^2 + dz^2];
    coef = prelel[[4]] * long / 2.;
    fe = Table[0.,{i,1,6}];
    fe[[3]] = coef;
    fe[[6]] = coef;
    Return[fe] ];
    
(*forces nodales dues à la déformation, qui fournissent le résidu d'équilibre par assemblage*)
Clear[res2];
res2[prelel_,xnel_,unel_]:=Module[{},

    dx = xnel[[2,1]] - xnel[[1,1]];
    dy = xnel[[2,2]] - xnel[[1,2]];
    dz = xnel[[2,3]] - xnel[[1,3]];
    long = Sqrt[dx^2 + dy^2 + dz^2];
    transf = {dx, dy, dz} / long;
    dl1 = transf[[1]]*unel[[1]] + transf[[2]]*unel[[2]] +
          transf[[3]]*unel[[3]];
    dl2 = transf[[1]]*unel[[4]] + transf[[2]]*unel[[5]] +
          transf[[3]]*unel[[6]];
    coef = prelel[[1]] * prelel[[2]] * (dl1 - dl2) / long;
    re = Table[0.,{i,1,6}];
    Do[ re[[i]] = transf[[i]] * coef; re[[i+3]] = - re[[i]],
      {i,1,3}];
    (* prise en compte des sollicitations réparties (seulement dans la direction "z" : poids propre ) *)
           If[prelel[[4]] != 0.0, 
         cc = prelel[[4]] * long / 2.;
         re[[3]] = re[[3]] - cc;
         re[[6]] = re[[6]] - cc, ];
    (* énergie de déformation *)
    energ = coef * (dl1 - dl2) / 2.;
    (*retour de l'energie de déformation et des forces nodales *)
    Return[{energ,re}] ];
    
(* calcul et impression des variables de l'élément liées aux gradients : 
                   déformation, effort et contrainte axiaux *)
Clear[grad2];
grad2[prelel_,xnel_,unel_]:=Module[{},

    dx = xnel[[2,1]] - xnel[[1,1]];
    dy = xnel[[2,2]] - xnel[[1,2]];
    dz = xnel[[2,3]] - xnel[[1,3]];
    long = Sqrt[dx^2 + dy^2 + dz^2];
    transf = {dx, dy, dz} / long;
    dl1 = transf[[1]]*unel[[1]] + transf[[2]]*unel[[2]] +
          transf[[3]]*unel[[3]];
    dl2 = transf[[1]]*unel[[4]] + transf[[2]]*unel[[5]] +
          transf[[3]]*unel[[6]];
    (* déformation axiale *)
           eps = (dl2 - dl1) / long;
    (* forces axiale *)
    fax = prelel[[1]] * prelel[[2]] * (dl2 - dl1) / long;
    (* contrainte *)
    sigma = fax / prelel[[2]];
    (* impression *)
    Print[" élément : ",iel];
    Print["      déformation  = ",eps];
    Print["      effort axial = ",fax];
    Print["      contrainte   = ",sigma];
    Return[{eps,fax,sigma}]  ];

(* tracé de l'élément avec déformée éventuelle *)
Clear[tracem2];
tracem2[xnel_,unel_,plot_]:=Module[{},

    (* coordonnées des nœuds de l'élément *)
    x1 = xnel[[1]]; x2 = xnel[[2]];
    (* tracé en "nplot"sous-éléments *)
    plot1=plot;
    plot1=Append[plot1,{Hue[0.3,1,1],Thickness[0.01],
          Line[{x1,x2}]} ];  
     Return[plot1] ];
    
(* tracé d'une grandeur définie aux noeuds par interpolation linéaire *)
Clear[tracei2];
tracei2[xnel_,unel_,ampli_,umin_,umax_,ncoupe_,uplot_,
        plot_]:=Module[{},

    (* coordonnées des nœuds de l'élément *)
    x1 = xnel[[1]] + ampli * unel[[1]];
    x2 = xnel[[2]] + ampli * unel[[2]];
    dx = (x2 - x1)/ncoupe;
    (* valeurs aux nœuds de l'élément de la fonction à tracer  *)
    upl1 = uplot[[1]]; upl2 = uplot[[2]];
    dupl = (upl2 - upl1)/ncoupe;
    
    
    (* tracé en "nplot"sous-éléments *)
    plot1=plot;
    x = x1; upl = upl1 + dupl/2.;
    hmax = 0.9; (* coefficient de réglage de la couleur, entre 0 et 1 *)
    Do[ 
      If[ umax!=umin,
           h = 1. - hmax + hmax *(upl - umin) / (umax - umin),
           h = 1.];
      plot1=Append[plot1,{Hue[h,1,1],Thickness[0.01],
                              Line[{x,x + dx}]} ];  
     x = x + dx; upl = upl + dupl,
      {i,1,ncoupe}];
    (* tracé *) 
    Return[plot1] ];

Elément de type 3 : triangle de Laplace linéaire à 2 dimensions et à 3 noeuds

(* propriétés élémentaires 
       prelel[[1]] : conductivité thermique
       prelel[[2]] : épaisseur
       prelel[[3]] : masse volumique
       prelel[[4]] : source par unité de surface *)
(* matrice de rigidité élémentaire *)
Clear[rgdt3];
rgdt3[prelel_,xnel_]:=Module[{},

    (* coordonnées des noeuds *)
    x1 = xnel[[1,1 ]]; x2 = xnel[[2,1 ]]; x3 = xnel[[3,1 ]];
    y1 = xnel[[1,2 ]]; y2 = xnel[[2,2 ]]; y3 = xnel[[3,2 ]];
    x21 = x2 - x1; x31 = x3 - x1; 
    y21 = y2 - y1; y31 = y3 - y1;
    (* centre de gravité *)
           xg = (x1 + x2 + x3)/3.;
    yg = (y1 + y2 + y3)/3.;
    aire2 = Abs[x21 * y31 - x31 * y21]; aire = aire2/2.;
    (* matrice jacobienne inverse *)
          jinv = {{y31,-y21},{-x31,x21}} / aire2;
    (* matrice B(2,3) *)
           matb = Table[0.,{i,1,2},{j,1,3}];
    ih = 2;  (* dimension de la matrice D *)
    Do[   matb[[i,1]] = -(jinv[[i,1]] + jinv[[i,2]]);
                   matb[[i,2]] = jinv[[i,1]];
          matb[[i,3]] = jinv[[i,2]],
      {i,1,ih}];
    (* matrice D *)
          coef = prelel[[1]] * prelel[[2]];
    matd = {{coef,0.},{0.,coef}};
    ke = Transpose[matb] . matd . matb * aire;
    Return[ke]  ];
    
    (* test de calcul d'une matrice rigidité élémentaire 
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    ke=rgdt3[prelel,xnel];
    Print[MatrixForm[ke]];
    Print[Sum[ke[[i,j]],{i,1,3},{j,1,3}]]; *)
    
(* matrice de rigidité tangente élémentaire *)
Clear[rtgt3];
rtgt3[prelel_,xnel_]:=Module[{},
   (* non implantée *)
    Return[ke] ];
    
(* matrice de masse élémentaire *)
Clear[mass3];
mass3[prelel_,xnel_]:=Module[{},

     (* coordonnées des noeuds *)
    x1 = xnel[[1,1 ]]; x2 = xnel[[2,1 ]]; x3 = xnel[[3,1 ]];
    y1 = xnel[[1,2 ]]; y2 = xnel[[2,2 ]]; y3 = xnel[[3,2 ]];
    x21 = x2 - x1; x31 = x3 - x1; 
    y21 = y2 - y1; y31 = y3 - y1;
    (* centre de gravité *)
           xg = (x1 + x2 + x3)/3.;
    yg = (y1 + y2 + y3)/3.;
    aire2 = Abs[x21 * y31 - x31 * y21]; aire = aire2/2.;
    cc = aire * prelel[[3]] / 12.;
    me = {{2.*cc,cc,cc},{cc,2.*cc,cc},{cc,cc,2.*cc}};
    Return[me] ];
    
    (* test de calcul d'une matrice masse
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    me=mass3[prelel,xnel];
    Print[MatrixForm[me]];
    Print[Sum[me[[i,j]],{i,1,3},{j,1,3}]]; *)
    
(* forces nodales dues à des sollicitations réparties *)
Clear[solr3];
solr3[prelel_,xnel_]:=Module[{},

    (* coordonnées des noeuds *)
    x1 = xnel[[1,1 ]]; x2 = xnel[[2,1 ]]; x3 = xnel[[3,1 ]];
    y1 = xnel[[1,2 ]]; y2 = xnel[[2,2 ]]; y3 = xnel[[3,2 ]];
    x21 = x2 - x1; x31 = x3 - x1; 
    y21 = y2 - y1; y31 = y3 - y1;
    (* centre de gravité *)
           xg = (x1 + x2 + x3)/3.;
    yg = (y1 + y2 + y3)/3.;
    aire2 = Abs[x21 * y31 - x31 * y21]; aire = aire2/2.;
    cc = prelel[[4]] * aire /3.;
    fe = {cc,cc,cc};
    Return[fe] ];
    
    (* test de calcul d'un vecteur de sollicitations élémentaire
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    fe=solr3[prelel,xnel];
    Print[MatrixForm[fe]];
    Print[Sum[fe[[i]],{i,1,3}]]; *)
    
(*forces nodales dues à la déformation, qui fournissent le résidu d'équilibre par assemblage*)
Clear[res3];
res3[prelel_,xnel_,unel_]:=Module[{},

    (* coordonnées des noeuds *)
    x1 = xnel[[1,1 ]]; x2 = xnel[[2,1 ]]; x3 = xnel[[3,1 ]];
    y1 = xnel[[1,2 ]]; y2 = xnel[[2,2 ]]; y3 = xnel[[3,2 ]];
    x21 = x2 - x1; x31 = x3 - x1; 
    y21 = y2 - y1; y31 = y3 - y1;
    (* centre de gravité *)
           xg = (x1 + x2 + x3)/3.;
    yg = (y1 + y2 + y3)/3.;
    aire2 = Abs[x21 * y31 - x31 * y21]; aire = aire2/2.;
    (* matrice jacobienne inverse *)
          jinv = {{y31,-y21},{-x31,x21}} / aire2;
    (* matrice B(2,3) *)
           matb = Table[0.,{i,1,2},{j,1,3}];
    ih = 2;
    Do[   matb[[i,1]] = -(jinv[[i,1]] + jinv[[i,2]]);
                   matb[[i,2]] = jinv[[i,1]];
          matb[[i,3]] = jinv[[i,2]],
       {i,1,ih}];
    (* matrice D *)
          coef = prelel[[1]] * prelel[[2]];
    matd = {{coef,0.},{0.,coef}};
    (* gradient *)
    grad = matb . unel;
    (* débit *)
    debit = matd . grad;
    (* énergie *)
    energ = 0.5 * (grad[[1]] * debit[[1]] + 
                   grad[[2]] * debit[[2]]) * aire;
    (* résidu élémentaire *)
    If[prelel[[4]] != 0.,
        coef = prelel[[4]] * aire / 3.,
        coef = 0.];
    re = Table[(matb[[1,i]] * debit[[1]] + 
                matb[[2,i]] * debit[[2]]) * aire - coef,{i,1,3}];            
                   
    Return[{energ,re}] ];
    
    (* test de calcul des forces nodales dues à la déformation 
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    unel  ={1.,1.,1.};
    result=res3[prelel,xnel,unel];
    Print[result];  *)
    
(* calcul et impression des variables de l'élément liées aux gradients : 
                   déformation, effort et contrainte axiaux *)
Clear[grad3];
grad3[prelel_,xnel_,unel_]:=Module[{},

    (* coordonnées des noeuds *)
    x1 = xnel[[1,1 ]]; x2 = xnel[[2,1 ]]; x3 = xnel[[3,1 ]];
    y1 = xnel[[1,2 ]]; y2 = xnel[[2,2 ]]; y3 = xnel[[3,2 ]];
    x21 = x2 - x1; x31 = x3 - x1; 
    y21 = y2 - y1; y31 = y3 - y1;
    (* centre de gravité *)
           xg = (x1 + x2 + x3)/3.;
    yg = (y1 + y2 + y3)/3.;
    aire2 = Abs[x21 * y31 - x31 * y21]; aire = aire2/2.;
    (* matrice jacobienne inverse *)
          jinv = {{y31,-y21},{-x31,x21}} / aire2;
    (* matrice B(2,3) *)
           matb = Table[0.,{i,1,2},{j,1,3}];
    ih = 2;
    Do[   matb[[i,1]] = -(jinv[[i,1]] + jinv[[i,2]]);
                   matb[[i,2]] = jinv[[i,1]];
          matb[[i,3]] = jinv[[i,2]],
       {i,1,ih}];
    (* matrice D *)
          coef = prelel[[1]] * prelel[[2]];
    matd = {{coef,0.},{0.,coef}};
    (* gradient *)
    grad = matb . unel;
    (* flux *)
    flux = matd . grad;
    norm = flux[[1]]^2 + flux[[2]]^2;
    (* impression *)
    Print[" élément : ",iel];
    Print["      gradient        : ",grad];
    Print["      flux en x et y  : ",flux];
    Print["      norme du flux   : ",norm];
    Return[{grad,flux,norm}]  ];
    
    (* test de calcul des variables liées au gradient  
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    unel  ={1.,1.,1.};
    result=grad3[prelel,xnel,unel]; *)
  
(* tracé de l'élément *)
Clear[tracem3];
tracem3[xnel_,unel_,plot_]:=Module[{},

    (* coordonnées des nœuds de l'élément *)
    x1 = xnel[[1]]; x2 = xnel[[2]]; x3 = xnel[[3]];
    (* tracé des 3 côtés de l'élément *)
    plot1=plot;
    plot1=Append[plot1,{Hue[0.3,1,1],Thickness[0.01],
                        Line[{x1,x2}]} ];
    plot1=Append[plot1,{Hue[0.3,1,1],Thickness[0.01],
                        Line[{x2,x3}]} ];
    plot1=Append[plot1,{Hue[0.3,1,1],Thickness[0.01],
                        Line[{x3,x1}]} ];
    Return[plot1] ];
     
    (* test de tracé d'un élément  
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    unel  ={1.,1.,1.};
    plot = {};
    plot = traceel3[xnel,unel,0.,plot];
    Show[Graphics[plot],Axes->True,AspectRatio->1.]; *)
    
(* tracé d'une grandeur définie aux noeuds par interpolation linéaire *)
Clear[tracei3];
tracei3[xnel_,unel_,ampli_,umin_,umax_,ncoupe_,uplot_,
        plot_,tr_,ntriangles_]:=Module[{},
    (* coordonnées des nœuds de l'élément *)
    x1 = xnel[[1]]; x2 = xnel[[2]]; x3 = xnel[[3]];
    (* valeurs aux nœuds de l'élément *)
    u1 = uplot[[1]];  u2 = uplot[[2]]; u3 = uplot[[3]];
     (* coefficient de réglage de la couleur, entre 0 et 1 *)
           hmax = 0.9; 
    plot1=plot;
    (* boucle sur les sous triangles *)
    Do[
      xx1=tr[[i,1,1]]*x1 + tr[[i,1,2]]*x2 + tr[[i,1,3]]*x3;
      xx2=tr[[i,2,1]]*x1 + tr[[i,2,2]]*x2 + tr[[i,2,3]]*x3;
      xx3=tr[[i,3,1]]*x1 + tr[[i,3,2]]*x2 + tr[[i,3,3]]*x3;
      (* valeur au centre du sous-triangle*)
      upl=tr[[i,4,1]]*u1 + tr[[i,4,2]]*u2 + tr[[i,4,3]]*u3;
      (* tracé du sous-élément *)
           If[ umax!=umin,
         h = 1. - hmax + hmax *(upl - umin) / (umax - umin),
         h = 1.];
    plot1=Append[plot1,{Hue[h,1.,1.], Polygon[{xx1,xx2,xx3}]} ]
   
    ,{i,1,ntriangles}]; 
    Return[plot1] ];
    
(* découpage d'un triangle de référence en sous-triangles 
                 (n par coté soit n couches de 1,2,3,..,n sous triangles) *)
Clear[soustriangles]
soustriangles[ncouches_] := Module[{},
    tr = {};
    imin=0.;jmin=0.;
    di=dj=1./ncouches;
    di3=di/3.;dj3=dj/3.;dj23=2.*dj/3.;

    Do[     (* boucle sur les couches *)
      i0=imin;
      j0=jmin;
      Do[   (* boucles sur les sous-triangles d'une couche *)
      If[ist==icouche,
       tr=Append[tr,{{i0,    j0,    1.-i0    -j0    },
                     {i0+di, j0,    1.-i0-di -j0    },
                     {i0,    j0+dj, 1.-i0    -j0-dj },
                     {i0+di3,j0+dj3,1.-i0-di3-j0-dj3}}],
                     
       tr=Append[tr,{{i0,    j0,    1.-i0    -j0},
                     {i0+di, j0,    1.-i0-di -j0},
                     {i0,    j0+dj, 1.-i0    -j0-dj},
                     {i0+di3,j0+dj3,1.-i0-di3-j0-dj3}}];
                      
       tr=Append[tr,{{i0,    j0,     1.-i0    -j0},
                     {i0,    j0+dj,  1.-i0    -j0-dj},
                     {i0-di, j0+dj,  1.-i0+di -j0-dj},
                     {i0-di3,j0+dj23,1.-i0+di3-j0-dj23}}]
               ];
       i0=i0-di;
       j0=j0+dj;
      ,{ist,1,icouche}];
     imin=imin+di,
    {icouche,1,ncouches}];
    ntriangles = Dimensions[tr][[1]];
    Return[{tr,ntriangles}] ];



     
    (* test de tracé d'un élément 
          ncoupe = 10; (* nombre de couches de sous-éléments par élément *)
           
    prelel={1.,1.,1.,1.};xnel={{0.,0.},{1.,0.},{0.,1.}};
    unel  ={1.,1.,1.}; 
    res=soustriangles[ncoupe];
    tr=res[[1]];ntriangles=res[[2]];
    uplot = {1.,0.,3.};
    umin = Min[uplot]; umax = Max[uplot];
    plot = {};
    plot = tracen3[xnel,unel,0.,umin,umax,ncoupe,uplot,plot,
                    tr,ntriangles];
   
    Show[Graphics[plot],
                Axes->True,AspectRatio->1.];  *)

Construction de la matrice globale de rigidité
   (penser à charger les fonctions élémentaires nécessaires)

(* création et initialisation à zéro de la matrice globale de rigidité *)
kg = Table[0.,{i,1,neq},{j,1,neq}];

(* boucle sur les éléments *)
Do[

(* extraction des propriétés de l'élément *)
prelel=prel[[1]]; 
(* à modifier pour prel[[iel]] si les propriétés ne sont pas constantes d'un élément à l'autre*)

(* liste des noeuds de l'élément *)
connel = conn[[iel]];
inel = Dimensions[connel][[1]]; (* nombre de noeuds de l'élément *)
idln = inel * ndln;             (* nombre de ddl de l'élément *)

(* construction de la localisation élémentaire *)
loce = Table[(connel[[in]]-1)*ndln+idl,{in,1,inel},{idl,1,ndln}];
loce = Flatten[loce]; (* table de localisation à 1 dimension *)

(* construction de la matrice de rigidité élémentaire *)
itpe = ntpe;    (* type d'élément *)
ke = Switch[itpe,
       1,rgdt1[prelel,xn[[connel]] ],
       2,rgdt2[prelel,xn[[connel]] ],
       3,rgdt3[prelel,xn[[connel]] ]  ];
If[print > 2,
Print[" élément : ",iel,", matrice de rigidité élémentaire :"];
Print[MatrixForm[ke]],];

(* assemblage de la matrice élémentaire dans la matrice globale *)
Do[             (* boucle sur les lignes   *)
  Do[           (* boucle sur les colonnes *)
kg[[loce[[i]],loce[[j]]]] = kg[[loce[[i]],loce[[j]]]] + ke[[i,j]],
     {j,1,idln}],
   {i,1,idln}],

{iel,1,nel}];   (* fin de la boucle d'assemblage *)

(* impression de la matrice globale *)
If[print > 2,
Print["matrice globale sans conditions aux limites:"];
Print[MatrixForm[kg]],]; 
(* impression de la diagonale de la matrice globale *)
If[print > 1,
diagkg = Table[0.,{i,1,neq}];
Do[diagkg[[i]]=kg[[i,i]],{i,1,neq}];
Print["diagonale matrice globale sans conditions aux limites:"];
Print[diagkg],];

Résultats exemple 1

Résultats exemple 3

Construction de la matrice globale de masse (si problème non stationnaire)
   (penser à charger les fonctions élémentaires)

(* création et initialisation à zéro de la matrice globale de masse *)
mg = Table[0.,{i,1,neq},{j,1,neq}];

(* boucle sur les éléments *)
Do[

(* extraction des propriétés de l'élément *)
prelel=prel[[1]]; 
(* à modifier pour prel[[iel]] si les propriétés ne sont pas constantes d'un élément à l'autre*)

(* liste des noeuds de l'élément *)
connel = conn[[iel]];
inel = Dimensions[connel][[1]]; (* nombre de noeuds de l'élément *)
idln = inel * ndln;             (* nombre de ddl de l'élément *)

(* construction de la localisation élémentaire *)
loce = Table[(connel[[in]]-1)*ndln+idl,{in,1,inel},{idl,1,ndln}];
loce = Flatten[loce]; (* table de localisation à 1 dimension *)

(* construction de la matrice de élémentaire *)
itpe = ntpe;    (* type d'élément *)
me = Switch[itpe,
       1,mass1[prelel,xn[[connel]] ],
       2,mass2[prelel,xn[[connel]] ],
       3,mass3[prelel,xn[[connel]] ]  ];
If[print > 2,
Print["--> élément : ",iel,", matrice de masse élémentaire :"];
Print[MatrixForm[me]],];

(* assemblage de la matrice élémentaire dans la matrice globale *)
Do[             (* boucle sur les lignes   *)
  Do[           (* boucle sur les colonnes *)
mg[[loce[[i]],loce[[j]]]] = mg[[loce[[i]],loce[[j]]]] + me[[i,j]],
     {j,1,idln}],
   {i,1,idln}],

{iel,1,nel}];   (* fin de la boucle d'assemblage *)

(* impression de la matrice globale *)
If[print > 2,
Print["matrice globale de masse:"];
Print[MatrixForm[mg]];
somme = Sum[mg[[i,j]],{i,1,neq},{j,1,neq}];
Print["somme des termes de la matrice de masse :  ",somme],];

Résultats exemple 1

Construction du vecteur global des sollicitations concentrées fg (second membre)

(* création et initialisation à zéro du vecteur  global des sollicitations concentrées *)
fgc = Table[0.,{i,1,neq}];

(* boucle sur les ddl sollicités *)
nsol = Dimensions[sollicit][[1]]; (* nombre de sollicitations *)
Do[
  in   = sollicit[[isol,1]];      (* numéro du noeud *)
  iddl = sollicit[[isol,2]];      (* numéro du ddl *)
  val  = sollicit[[isol,3]];      (* valeur de la sollicitation *)
   
  ieq  = (in - 1) * ndln + iddl;  (* numéro global du ddl sollicité *)
  fgc[[ieq]] = fgc[[ieq]] + val;

,{isol,1,nsol}];

(* création et initialisation à zéro du vecteur  global des sollicitations (second membre) *)
fg = fgc;

(* impression du vecteur second membre *)
If[print > 1,
Print["vecteur second membre sans conditions aux limites:"];
Print[fg],];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Construction du vecteur global des sollicitations réparties fg (second membre)

(* initialisation du vecteur global des sollicitations, s'il n'existe pas déjà *)
fgr = Table[0.,{i,1,neq}];
(* début de la boucle sur les éléments  *)
Do[
(* liste des noeuds de l'élément *)
connel = conn[[iel]];

(* extraction des propriétés de l'élément *)
prelel=prel[[1]]; (* à modifier pour Prel[[iel]]si les propriétés ne 
                        sont pas constantes d'un élément à l'autre*)
(* calcul des sollicitations concentrées élémentaires *)
fe = Switch[itpe,
       1,solr1[prelel,xn[[connel]]],
       2,solr2[prelel,xn[[connel]]],
       3,solr3[prelel,xn[[connel]]]  ];
If[print > 1,
Print["élément : ",iel,
      ", sollicitations réparties élémentaire : "];
Print[fe],];

(* construction de la localisation élémentaire *)
inel = Dimensions[connel][[1]]; (* nombre de noeuds de l'élément *)
idln = inel * ndln;             (* nombre de ddl de l'élément *)
loce = Table[(connel[[in]]-1)*ndln+idl,{in,1,inel},{idl,1,ndln}];
loce = Flatten[loce];           (* vecteurde localisation à 1 dimension *)

(* assemblage du résidu élémentaire dans le résidu global *)
    Do[             (* boucle sur les lignes   *)
       fgr[[loce[[i]]]] = fgr[[loce[[i]]]] + fe[[i]],
      {i,1,idln}],

 {iel,1,nel}];(* fin de la boucle sur les éléments *)
 
 fg = fg + fgr;
 
(* impression du vecteur résidu global *)
If[print > 1,
Print[" sollicitations réparties : "];
Print[fgr];
Print["somme des sollicitations réparties"];
Print[Sum[fgr[[i]],{i,1,neq}]];
Print["sollicitations totales : "];
Print[fg],];

Résultats exemple 3

Calcul des valeurs propres et vecteurs propres de  la matrice globale avant application des conditions aux limites (facultatif)

(* utiliser la méthode de calcul des valeurs et vecteurs propres de Mathematica *)
eig = Eigensystem[kg];

(* impression *)
If[print > 1,
Print["   Avant prise en compte des conditions aux limites :"];
Print["valeurs propres : ", eig[[1]]],];
If[print > 2,
Print["vecteurs propres : "];
Print[MatrixForm[eig[[2]]]],];

Application des conditions aux limites

(* boucle sur les ddl sollicités *)
ncond = Dimensions[condlim][[1]]; (* nombre de conditions aux limites *)
Do[
  in   = condlim[[icond,1]];      (* numéro du noeud *)
  iddl = condlim[[icond,2]];      (* numéro du ddl   *)
  val  = condlim[[icond,3]];      (* valeur imposée *)
   
  ieq  = (in - 1) * ndln + iddl;  (* numéro global du ddl imposé *)
  
  Do[
      (* retrancher au second membre la colonne de la matrice globale  multipliée par la valeur imposée *)
      fg[[i]] = fg[[i]] - val * kg[[i,ieq]];
      (* annuler la ligne et la colonne de la matrice globale *)
      kg[[i,ieq]] = 0.;
      kg[[ieq,i]] = 0.,
  {i,1,neq}];
  
  (* mettre un 1 sur la diagonale de la matrice globale *)
  kg[[ieq,ieq]] = 1.;
  (* mettre la valeur imposée en ligne "icond" du vecteur second membre *)
  fg[[ieq]]       = val,

{icond,1,ncond}];

(* impression de la matrice globale *)
If[print > 2,
Print["matrice globale avec conditions aux limites:"];
Print[MatrixForm[kg]],];

(* impression du vecteur second membre *)
If[print > 1,
Print["vecteur second membre avec conditions aux limites:"];
Print[fg],];

Résultats exemple 1

Résultats exemple 3

Calcul des valeurs propres et vecteurs propres de  la matrice globale après application des conditions aux limites (facultatif)

(* utiliser la méthode de calcul des valeurs et vecteurs propres de Mathematica *)
eig = Eigensystem[kg];

(* impression *)
If[print > 1,
Print["   Aprés prise en compte des conditions aux limites :"];
Print["valeurs propres : ", eig[[1]]],];
If[print > 2,
Print["vecteurs propres : "];
Print[MatrixForm[eig[[2]]]],];

Résultats exemple 1

Etape 3 : résolution

Résolution du système d'équations linéaires : kg . u = fg

(* utilisation de la méthode de résolution directe de Mathematica *)
un = LinearSolve[kg,fg];

(* impressions *)
If[print > 2,
Print["matrice globale : "];
Print[MatrixForm[kg]],];
If[print > 1,
Print["vecteur second membre : "];
Print[fg],];
If[print > 0,
Print["solution : "];
Print[un],];

Résultats exemple 1

Résultats exemple 3

Résolution du problème de valeurs et vecteurs propres :  
              (kg - λ * mg) . u = 0

(* utiliser la méthode de calcul des valeurs et vecteurs propres de Mathematica *)
eig = Eigensystem[Inverse[mg] . kg];

(* impression *)
If[print > 1,
Print["  Problème de valeurs et vecteurs propres"];
Print["        (kg - λ * mg) . u = 0 :"];
Print["valeurs propres : ", eig[[1]]],];
If[print > 2,
Print["vecteurs propres : "];
Print[MatrixForm[eig[[2]]]],];

Résolution d'un problème non linéaire par la méthode de Newton Raphson

(* entreé du nombre de pas *)
npas = 3;
(*entrée de la taille de chaque pas *)
dcharge = 1.;
(*entrée du nombre maximum d'itérations par pas *)
niter = 2;
(*entrée des tolérance sur les inconnues et sur le résidu *)
epsdun = 10.^-3;
epsres = 10.^-2;

(* initialiser les inconnues, ici à zéro *)
un = Table[0.,{i,1,neq}];
(* initialiser à zéro le niveau de chargement *)
charge = 0.;
(* boucle sur les pas *)
Do[
(* initialiser les inconnues, ici à zéro *)
charge = charge + dcharge;
Print["==============================   "];
Print["niveau de chargement : ",charge];
Print["==============================   "];
Print[" "];

   (* boucle sur les itérations *)
   Do[
   Print["------------------"];
   Print["itération : ",iter];
   Print["----------------- "];
   (* assemblage de la matrice tangente et du vecteur des résidus *)
             (* initialisation du vecteur global des résidus et de la matrice dangente *)
     rg = Table[0.,{i,1,neq}];
     kt = Table[0.,{i,1,neq},{j,1,neq}];
     (* début de la boucled'assemblage sur les éléments  *)
     Do[
     (* liste des noeuds de l'élément *)
     connel = conn[[iel]];
     (* extraction des propriétés de l'élément *)
     prelel=prel[[1]]; (* à modifier pour Prel[[iel]]si les 
           propriétés ne sont pas constantes d'un élément à l'autre*)
     (* extraction des ddl d'un élément d'un élément *)
              unel = Flatten[Table[un[[jd + (connel[[in]] - 1) *
                    ndln]], {in,1,nnel},{jd,1,ndln}]];
     (* calcul du résidu et de l'énergie élémentaires *)
             result = Switch[itpe,
       1,res1[prelel,xn[[connel]],unel],
       2,res2[prelel,xn[[connel]],unel],
       3,res3[prelel,xn[[connel]],unel]  ];
     re = result[[2]];
     (* calcul de la matrice tangente
             kte = Switch[itpe,
       1,rgtg1[prelel,xn[[connel]],unel],
       2,rgtg2[prelel,xn[[connel]],unel],
       3,rgtg3[prelel,xn[[connel]],unel]  ]; *)
       
             (* calcul de la matrice rigidité (Méthode de Newton modifiée ) *)
             kte = Switch[itpe,
       1,rgdt1[prelel,xn[[connel]]],
       2,rgdt2[prelel,xn[[connel]]],
       3,rgdt3[prelel,xn[[connel]]]  ];

     (* construction de la localisation élémentaire *)
             inel = Dimensions[connel][[1]]; (* nombre de noeuds de l'élément *)
     idln = inel * ndln;             (* nombre de ddl de l'élément *)
     loce = Table[(connel[[in]]-1)*ndln+idl,{in,1,inel},
                                           {idl,1,ndln}];
     loce = Flatten[loce];           (* vecteurde localisation à 1 dimension *)

     (* assemblage du résidu élémentaire dans le résidu global *)
     Do[             (* boucle sur les lignes   *)
          rg[[loce[[i]]]] = rg[[loce[[i]]]] + re[[i]];
          Do[          (* boucle sur les colonnes *)
          kt[[loce[[i]],loce[[j]]]] = kt[[loce[[i]],loce[[j]]]]+ 
               kte[[i,j]],
            {j,1,idln}],
       {i,1,idln}],

    {iel,1,nel}];   (* fin de la boucle sur les éléments *)
            (* prise en compte des sollicitations au niveau de chargement atteint *)
            rg = charge * fg - rg;
            (* impression du vecteur résidu global de deséquilibre *)
            If[print > 1,
     Print["   (avant  conditions aux limites)"];
     Print[" "];
             Print[" résidu global de déséquilibre et réactions : "];
     Print[rg],];
     (* impression de la matrice tangente globale *)
             If[print > 2,
     Print[" matrice tangente globale : "];
     Print[kt],];
     
     (* application des conditions aux limites à la matrice tangente et au résidu global *)
       (* boucle sur les ddl sollicités *)
     ncond = Dimensions[condlim][[1]]; (* nombre de conditions 
                                                                                                             aux limites *)
     Do[
       in   = condlim[[icond,1]];      (* numéro du noeud *)
       iddl = condlim[[icond,2]];      (* numéro du ddl   *)
       val  = condlim[[icond,3]];      (* valeur imposée au ddl, de référence *)
   
       ieq  = (in - 1) * ndln + iddl;  (* numéro global du ddl imposé *)
                   (* valeur qu'il faut imposer à l'accroissement de ddl 
                             (condition à la limite supposée ici évoluer proportionnellement au paramètre de charge) 
                             ... on peut faire une autre hypothèse ...  *)
                   valdu = charge * val - un[[ieq]];
  
       Do[
          (* retrancher au second membre la colonne de la matrice globale  
                                                                                            multipliée par la valeur imposée *)
          rg[[i]] = rg[[i]] - valdu * kt[[i,ieq]];
          (* annuler la ligne et la colonne de la matrice globale *)
          kt[[i,ieq]] = 0.;
          kt[[ieq,i]] = 0.,
         {i,1,neq}];
       (* mettre un 1 sur la diagonale de la matrice globale *)
       kt[[ieq,ieq]] = 1.;
       (* mettre la valeur imposée en ligne "icond" du vecteur second membre *)
       rg[[ieq]] = valdu,
       {icond,1,ncond}];
       
     If[print > 1,
     Print["   (après  conditions aux limites)"];
     Print[" "];
             Print[" résidu global de déséquilibre et réactions : "];
     Print[rg],];
     (* impression de la matrice tangente globale *)
             If[print > 2,
     Print[" matrice tangente globale : "];
     Print[kt],];
       
   (* résolution pour calculer l'accroissement des inconnues *)
       dun = LinearSolve[kt,rg];
       
   (* mise à jour des inconnues *)
   un = un + dun;
   
   If[print > 1,
      Print[" "];
      Print["accroissement des inconnues = "];
      Print[dun];
      Print["inconnues = "];
      Print[un],];
      Print[" "];
      
   (* analyse de la convergence : calcul des normes des accroissements de ddl et des résidu *)
   normdun = Sqrt[Sum[dun[[i]]^2,{i,1,neq}]];
   normres = Sqrt[Sum[rg[[i]]^2,{i,1,neq}]];
   If[print > 0,
   Print["norme des accroissements des inconnues = ",
          normdun," (eps_du  = ",epsdun,")"];
   Print["norme du  résidu global                = ",
          normres," (eps_res = ",epsres,")"],
     ];
     
   (* convergence *)
   If[normdun < epsdun && normres < epsres,converge = True;
      Print[" "];
      Print["---> CONVERGENCE EN ",iter," ITERATIONS"];
      Break[],converge = False],
      
   {iter,1,niter}];  (* fin de boucle sur les itérations de résolution *)
         (*non convergence *)
     If[converge == False,
     Print[" "];
     Print["--->NON CONVERGENCE EN ",niter," ITERATIONS"],];
     
If[print > 0,
Print[" "];
Print["inconnues = "];
Print[un],],
{ipas,1,npas}];            (* fin de boucle sur les pas de chargement *)

Résultats exemple 3

Intégration du système différentiel par rapport au temps du second ordre par
     la méthode de Newmark :
                   mg . u°° + cg . u° + kg . u = fg(t)
                 où  :  u°°(t) est l'accélération des nœuds
                         u° (t)  est la vitesse des nœuds
                         u  (t)  est la position des nœuds
    (on suppose les matrice globales de rigidité kg et de masse mg déjà calculées)

(* entrée de paramètres variés *)
a = 0.5; b = 0.5;         (* paramètres de la méthodes de Newmark *)
deltat = 0.5;             (* pas de temps *)
npas = 31;                (* nombre de pas d'intégration *)
lambda = 1.; mu = 0.;     (* coefficients de définition de la matrice d'amortissement *)

tmax = deltat * npas;     (* temps maximal *)
ntemps = npas + 1;        (* nombre d'instants étudiés incluant l'instant initial *)

(*définir l'évolution du vecteur second membre au cours du temps *)
Clear[fgt];
omega = N[Pi]/10.;
fgt[t_]:=fg * Sin[omega * t];  (* on peut compliquer à plaisir ... *) 
(* conditions initiales à l'instant t = 0, à modifier si nécessaire *)
u   = Table[0.,{i,1,neq}];    (* déplacements initiaux *)
up  = Table[0.,{i,1,neq}];    (* vitesse initiale *)

(* matrice d'amortissement définie ici comme une combinaison linéaire des matrices de masse et de rigidité *)
cg = lambda * kg + mu * mg;

(* calculer l'accélération initiale grâce aux équations d'équilibre à l'instant t = 0 *)
upp = LinearSolve[mg,fgt[0.] - (cg . up) - (kg . u)];

(* tableaux de stockage des déplacements et réactions  pour les tracés ultérieurs *)
uplot    = Table[0.,{i,1,ntemps},{j,1,neq}];
upplot   = Table[0.,{i,1,ntemps},{j,1,neq}];
uppplot  = Table[0.,{i,1,ntemps},{j,1,neq}];
plotreac = Table[0.,{i,1,ntemps},{j,1,neq}];
(* stocker la solution initiale *)
uplot[[1]]    = u;
upplot[[1]]   = up;
uppplot[[1]]  = upp;
plotreac[[1]] = mg . upp + cg . up + kg . u;
(* coefficients de calcul de k1 et r1  évalués à l'instant t+deltat *)
a1=deltat * a; b1=deltat*deltat*b/2.;
b2=deltat * deltat * (1. - b)/2.;
a2=deltat * a; a3=deltat * deltat * (2.* a - b)/2.; 
a4=deltat * deltat * deltat * (a-b)/2.;
(* coef de calcul de la solution à  l'instant t+deltat *)
b3=2./(b*deltat*deltat); b4=2./(b*deltat); b5=(1./b)-1.;
a5=1.-a; 
b6=deltat*deltat*(1.- b)/2; b7=deltat*deltat*b/2.;

(* boucle d'intégration *)
t=0.;
Do[

(* solution à l'instant t+deltat *)
k1 = mg + a1 * cg + b1 * kg;
r1 = b1 * fgt[t+deltat] + mg . (u + deltat * up + b2 * upp) +
     cg . (a2 * u + a3 * up + a4 * upp);
  (*   
(* application des conditions aux limites à k1 et r1 *)
(* boucle sur les ddl sollicités *)
ncond = Dimensions[condlim][[1]]; (* nombre de conditions aux limites *)
Do[
  in   = condlim[[icond,1]];      (* numéro du noeud *)
  iddl = condlim[[icond,2]];      (* numéro du ddl   *)
  val  = condlim[[icond,3]];      (* valeur imposée *)
   
  ieq  = (in - 1) * ndln + iddl;  (* numéro global du ddl imposé *)
  (* boucle sur les lignes *)
  Do[
      (* retrancher du second membre la colonne de la matrice globale  multipliée par la valeur imposée *)
      r1[[i]] = r1[[i]] - val * k1[[i,icond]];
      (* annuler la ligne et la colonne de la matrice globale *)
      k1[[i,icond]] = 0.;
      k1[[icond,i]] = 0.,
  {i,1,neq}];
  
  (* mettre un 1 sur la diagonale de la matrice globale *)
  k1[[icond,icond]] = 1.;
  (* mettre la valeur imposée en ligne "icond" du vecteur second membre *)
  r1[[icond]]       = val,

{icond,1,ncond}];
    *)
(* calcul des déplacements à l'instant t + deltat *)
u1 = LinearSolve[k1,r1];

(* calcul de l'accélération à l'instant t+deltat *)
u1pp = b3 * (u1 - u) - b4 * up - b5 * upp;

(* calcul de la vitesse à l'instant t + deltat *)
u1p  = up + deltat * (a5 * upp + a * u1pp);

(* mise à jour de la solution et sauvegarde *)
t = t + deltat;
u   = u1;
up  = u1p;     
upp = u1pp;

(* sauvegarde des deplacements pour le tracé *)
uplot[[ip+1]]   = u;
upplot[[ip+1]]  = up;
uppplot[[ip+1]] = upp;

(* calcul et sauvegares des réactions pour le tracé *)
reac = mg . upp + cg . up + kg . u;
plotreac[[ip]] = reac,

  {ip,1,npas}];    (* fin de boucle d'intégration *)
  
(* tracé de la variation du ddl numéro iddl au cours des "npas" pas *)
iddl = 1;
ListPlot[uplot[[Table[i,{i,1,ntemps}],iddl]],
      PlotRange->All,PlotJoined->True];
iddl = 2;
ListPlot[uplot[[Table[i,{i,1,ntemps}],iddl]],
      PlotRange->All,PlotJoined->True];
ListPlot3D[uplot,PlotRange->All];

Résultats exemple 1

Etape 4 : post-traitement

Impression de la solution aux noeuds

Print["solution : ",un];

Résultats exemple 1

Résultats exemple 3

Impression des gradients (inconnues secondaires) de chaque élément
   (penser à charger les fonctions élémentaires)

(* début de la boucle sur les éléments  *)
Do[
(* liste des noeuds de l'élément *)
connel = conn[[iel]];

(* extraction des propriétés de l'élément *)
prelel=prel[[1]]; (* à modifier pour Prel[[iel]]si les propriétés ne 
                        sont pas constantes d'un élément à l'autre*)
(* extraction des ddl d'un élément d'un élément *)
unel = Flatten[Table[un[[jd + (connel[[in]] - 1) * ndln]],
                     {in,1,nnel},{jd,1,ndln}]];
(* calcul des gradients d'un élément *)
result = Switch[itpe,
       1,grad1[prelel,xn[[connel]],unel],
       2,grad2[prelel,xn[[connel]],unel],
       3,grad3[prelel,xn[[connel]],unel]  ],

{iel,1,nel}];    (* fin de la boucle sur les éléments *)

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Assemblage du vecteur de résidu de déséquilibre incluant les réactions
   (penser à charger les fonctions élémentaires)

(* initialisation du vecteur global des résidus *)
rg = Table[0.,{i,1,neq}];
(* début de la boucle sur les éléments  *)
Do[
(* liste des noeuds de l'élément *)
connel = conn[[iel]];

(* extraction des propriétés de l'élément *)
prelel=prel[[1]]; (* à modifier pour Prel[[iel]]si les propriétés ne 
                        sont pas constantes d'un élément à l'autre*)
(* extraction des ddl d'un élément d'un élément *)
unel = Flatten[Table[un[[jd + (connel[[in]] - 1) * ndln]],
                     {in,1,nnel},{jd,1,ndln}]];
(* calcul du résidu et de l'énergie élémentaires *)
result = Switch[itpe,
       1,res1[prelel,xn[[connel]],unel],
       2,res2[prelel,xn[[connel]],unel],
       3,res3[prelel,xn[[connel]],unel]  ];
energ = result[[1]]; re = result[[2]];

(* construction de la localisation élémentaire *)
inel = Dimensions[connel][[1]]; (* nombre de noeuds de l'élément *)
idln = inel * ndln;             (* nombre de ddl de l'élément *)
loce = Table[(connel[[in]]-1)*ndln+idl,{in,1,inel},{idl,1,ndln}];
loce = Flatten[loce];           (* vecteurde localisation à 1 dimension *)

(* assemblage du résidu élémentaire dans le résidu global *)
     Do[             (* boucle sur les lignes   *)
          rg[[loce[[i]]]] = rg[[loce[[i]]]] + re[[i]],
       {i,1,idln}],

 {iel,1,nel}];   (* fin de la boucle sur les éléments *)
 
(* impression du vecteur résidu global de deséquilibre *)
Print[" résidu global de déséquilibre et réactions : "];
Print[rg - fgc];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Tracer du maillage
   (penser à charger les fonctions élémentaires)

(* initialisation de la table de tracer *)
plot = {};

(* début de la boucle de tracé virtuel des éléments  *)
Do[
(* liste des noeuds de l'élément *)
connel = conn[[iel]];
(* extraction des ddl d'un élément d'un élément *)
unel = Flatten[Table[un[[jd + (connel[[in]] - 1) * ndln]],
                     {in,1,nnel},{jd,1,ndln}]];
(* ajout d'un élément dans la table de tracer *)
plot = Switch[itpe,
       1,tracem1[xn[[connel]],unel,plot],
       2,tracem2[xn[[connel]],unel,plot],
       3,tracem3[xn[[connel]],unel,plot]  ],

{iel,1,nel}];   (* fin de la boucle de construction de la table de tracer *)

(* début de la boucle de tracé virtuel des nœuds  *)
Do[
(* table des ddl du nœuds *)
unn = Table[un[[i + (in - 1) * ndln]],{i,1,ndln}];
(* ajout d'un nœud dans la table de tracer *)
Switch[ndim,
    (* cas d'un problème à 1 dimension dans lequel les coordonnées "y" ne sont pas définies *)
    1,xxn = {xn[[in]][[1]],0.};
      plot = Append[plot,{Hue[1,1,1],PointSize[0.02],
                          Point[xxn]} ],
    (* cas d'un problème à 2 dimensions *)
    2,xxn = xn[[in]];
      plot = Append[plot,{Hue[1,1,1],PointSize[0.02],
                          Point[xxn]} ],
    (* cas d'un problème à 3 dimensions *)
    3,xxn = xn[[in]];
      plot = Append[plot,{Hue[1,1,1],PointSize[0.02],
                          Point[xxn]} ] ],
 {in,1,nnt}];   (* fin de la boucle de construction de la table de tracer *)

(* tracé effectif des éléments et de nœuds *)
Switch[ndim,
          1,Show[Graphics[plot],Axes->True],
          2,Show[Graphics[plot],Axes->True],
          3,Show[Graphics3D[plot],Axes->True] ];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3

Tracer des iso-valeurs d'une grandeur définie aux nœuds sur un maillage éventuellement déformé :
   (penser à charger les fonctions élémentaires)

(* coefficient d'amplification de la déformée (nul pour une simple tracé du maillage) *)
ampli = 0.;
(* nombre de sous-éléments utilisés pour le tracé *)
ncoupe = 10;
(* définition de la grandeur définie aux nœuds, à tracer sur les éléments *)
uplot = Table[un[[1 + (in-1)*ndln]],{in,1,nnt}];  (* ici déplacement "u" *)


(* initialisation de la table de tracer *)
plot = {};
(* valeurs minimum et maximum des valeurs à tracer *)
umin = Min[uplot]; umax = Max[uplot];

(* début de la boucle de tracé virtuel des éléments  *)
Do[
(* liste des noeuds de l'élément *)
connel = conn[[iel]];
(* extraction des ddl d'un élément d'un élément *)
unel = Table[un[[jd + (connel[[in]] - 1) * ndln]],
                     {in,1,nnel},{jd,1,ndln}];
(* ajout d'un élément dans la table de tracer *)
plot = Switch[itpe,
       1,tracei1[xn[[connel]],unel,ampli,umin,umax,ncoupe,
                 uplot[[connel]],plot],
       2,tracei2[xn[[connel]],unel,ampli,umin,umax,ncoupe,
                 uplot[[connel]],plot],
       3,If[iel == 1,res=soustriangles[ncoupe];
                     tr=res[[1]];ntriangles=res[[2]],];
         tracei3[xn[[connel]],unel,ampli,umin,umax,ncoupe,
                 uplot[[connel]],plot,tr,ntriangles]  ],

{iel,1,nel}];   (* fin de la boucle de construction de la table de tracer *)

(* début de la boucle de tracé virtuel des nœuds  *)
Do[
(* table des ddl du nœuds *)
unn = Table[un[[i + (in - 1) * ndln]],{i,1,ndln}];
(* ajout d'un nœud dans la table de tracer *)
Switch[ndim,
    (* cas d'un problème à 1 dimension dans lequel les coordonnées "y" ne sont pas définies *)
    1,If[ampli == 0.,xxn = {xn[[in]][[1]],0.},
              xxn = {xn[[in]][[1]] + ampli * un[[in]] ,0.}];
      plot = Append[plot,{Hue[1,1,1],PointSize[0.02],
                          Point[xxn]} ],
    (* cas d'un problème à 2 dimensions *)
    2,If[ampli == 0.,xxn = xn[[in]] ,
              xxn = xn[[in]] + ampli * unn ];
      plot = Append[plot,{Hue[1,1,1],PointSize[0.02],
                          Point[xxn]} ],
    (* cas d'un problème à 3 dimensions *)
    3,If[ampli == 0.,xxn = xn[[in]] ,
              xxn = xn[[in]] + ampli * unn ];
      plot = Append[plot,{Hue[1,1,1],PointSize[0.02],
                          Point[xxn]} ] ],

 {in,1,nnt}];   (* fin de la boucle de construction de la table de tracer *)


(* tracé effectif des éléments et de nœuds *)
Switch[ndim,
          1,Show[Graphics[plot],Axes->True],
          2,Show[Graphics[plot],Axes->True],
          3,Show[Graphics3D[plot],Axes->True] ];

Résultats exemple 1

Résultats exemple 2

Résultats exemple 3