************************************ calcul1.dgibi
******** optiON de calcul :
********** dimeSION 2 , elemENT TRI3 , mode de contRAINTE planE 
********** echo 0 (pour peu d'affichage) ou echo 1 (pour beaucoup d'affichage)
********** densITE du maillage exprimée dans l'unité des cotes
opti dime 2 elem tri3 mode plan cont echo 1 ;
********........ cotes de la structure
l0 = 10.0 ;
h0 = 5.0 ;
dens 1.50 ;
********........ géométrie
***** .......... les points définissant le contour extérieur
p1 = 0.0 0.0 ;
p2 = l0 0.0 ;
p3 = l0 h0 ;
p4 = 0.0 h0 ;
c1 = (l0/2.0) (h0/2.0) ;
*****............ les lignes définissant le contour extérieur
l1 = droi p1 p2 ;
l2 = droi p2 p3 ;
l3 = cerc p3 c1 p4 ;
l4 = droi p4 p1 ;
*****............ le contour extérieur
cont0 = l1 et l2 et l3 et l4 ;
*** on trace à l'écran
trac cont0 titr 'Le contour' ;
********......... maillage d'une surface plane
mail0 = surf cont0 plan ;
trac mail0 titr 'Maillage' ;
*****............ on compte le nombre de noeud
nbno0 = nbno mail0 ;
mess 'Nombre de noeuds :' nbno0 ;
*****............ on compte le nombre d'élément
nbel0 = nbel mail0 ;
mess 'Nombre d éléments :' nbel0 ;
********.......... modèle (homogène et linéaire sous entendu)
mod0 = mode mail0 mecanique elastique isotrope ;
**** caractéristiques du matériau (module de Young youn et coef. de Poisson nu)
mat0 = mate mod0 youn 210000.0 nu 0.30 ;
********.......... conditions limites : blocage du déplacement sur la ligne l1
cl0 = bloq depla l1 ;
*********......... chargement ponctuel
fp1 = -1000.0 ;
***............... suivant y au point p4
for1 = forc fy fp1 p4 ;
***............... amplitude de représentation à l'écran
amp0 = 1.0e-2 ;
***............... on créé un vecteur pour la représentation à l'écran
vect1 = vect for1 amp0 fx fy bleu ;
********.......... chargement réparti
fr1 = 1000.0/h0 ;
****.......... vecteur force linéique 
vfr = fr1 0.0 ;
**** ou force surfacique (car l'épaisseur de la structure est de 1 dans l'unité utilisée) imposée sur l2
for2 = fsur mass mod0 l2 vfr ;
***............... on créé un autre vecteur pour la représentation à l'écran
vect2 = vect for2 amp0 fx fy bleu ;
trac (vect2 et vect1) mail0 titr 'Chargements' ;
**** for0 contient toutes les actions extérieures
for0 = for1 et for2 ;
******** matrice de rigidité
rig0 = rigi mod0 mat0 ;
******** prise en compte des conditions limites
******** dans la matrice de raideur
rig0 = rig0 et cl0 ;
********************* résolution d'un système de nbno0*nbinc équations (où ici nbinc=2 (u et v))
******** [K][x]=[F] avec
******** [K] = rig0 ; [F] = for0 et [x] = res0 contient les ddl (u et v)
******** en chacun des noeuds du maillage
res0 = reso rig0 for0 ;
temps ;
********
def0 = defo res0 mail0 0.0 ;
******** déformée (amplitude calculée automatiquement)
def1 = defo res0 mail0 rouge ;
trac (def0 et def1) titre 'déformée' ;
******** déformée (amplitude vraie)
def1 = defo res0 mail0 1.0 rouge ;
trac (def0 et def1) titre 'déformée' ;
******** valeur déplacement en P4
mess '' ;
****............ 1ère et 2nde coordonnées du point p4
xa = coor 1 p4 ; ya = coor 2 p4 ;
mess 'x,y en P4 : 'xa ya ;
****............ extraction de la solution (res0) du déplacement suivant x (ux) [ou y (uy)] au point p4
ux4 = extr res0 ux p4 ;
mess 'Déplacement suivant x en P4 : 'ux4 ;
uy4 = extr res0 uy p4 ;
mess 'Déplacement suivant y en P4 : 'uy4 ; mess '' ;
******** valeur déplacement en (l0/2,h0)
****..........on cherche le point du maillage le plus proche de ces coordonnées et on l'appelle p5
p5 = mail0 poin proc (l0/2.0 h0) ;
mess 'x,y P5 : ' (l0/2.0) h0 ;
x5 = coor 1 p5 ; y5 = coor 2 p5 ;
mess 'x,y le plus proche de P5 : 'x5 y5 ;
ux5 = extr res0 ux p5 ;
mess 'Déplacement suivant x en P5 : 'ux5 ;
uy5 = extr res0 uy p5 ;
mess 'Déplacement suivant y en P5 : 'uy5 ; mess '' ;
******** action aux liaisons
****............. calculée
reac0 = reac rig0 res0 ;
****............. représentée
vreac0 = vect reac0 1.0e-2 fx fy bleu ;
trac vreac0 mail0 titr 'Action aux liaisons';
****............. résultante de ces actions 
fg0 = resu reac0 ;
mess 'Force globale :' ; list fg0 ;
**********************************************
******** calcul des contraintes dans l'élément
sig0 = sigm mod0 mat0 res0 ;
******** moyenne des contraintes aux noeuds
signo0 = chan chpo mod0 sig0 ;
*
******** tracé de sigmaxx sur les elements
sxx0 = exco smxx sig0 ;
trac mod0 sxx0 titr 'Sigma xx' ;
******** tracé de sigmaxx moyenné aux noeuds
sxxno0 = exco smxx signo0 ;
trac mail0 sxxno0 titr 'Sigma xx' ;
******** tracé de sigmaxy sur les elements
sxy0 = exco smxy sig0 ;
trac mod0 sxy0 titr 'Sigma xy' ;
******** tracé de sigmaxy moyenné aux noeuds
sxyno0 = exco smxy signo0 ;
trac mail0 sxyno0 titr 'Sigma xy' ;
******** tracé de sigmayy sur les elements
syy0 = exco smyy sig0 ;
trac mod0 syy0 titr 'Sigma yy' ;
******** trace de sigmayy moyenné aux noeuds
syyno0 = exco smyy signo0 ;
trac mail0 syyno0 titr 'Sigma yy' ;
***** avec les mêmes isovaleurs
isov0 = prog -550 -500 -450 -400 -350 -300
       -250 -200 -150 -100 -50 0 50 100 150 200 250 ;
trac mail0 sxxno0 titr 'Sigma xx' isov0  ;
trac mail0 sxyno0 titr 'Sigma xy' isov0  ;
trac mail0 syyno0 titr 'Sigma yy' isov0  ;
***********************************************
******** calcul des déformations dans l'élément
eps0 = epsi mod0 res0 ;
******** moyenne des déformations aux noeuds
epsno0 = chan chpo mod0 eps0 ;
******** tracé de epsilonxx moyenné aux noeuds
exxno0 = exco epxx epsno0 ;
trac mail0 exxno0 titr 'Epsilon xx' ;
******** tracé de epsilonxy moyenné aux noeuds
exyno0 = exco gaxy epsno0 ;
trac mail0 exyno0 titr 'Epsilon xy' ;
******** tracé de epsilonyy moyenné aux noeuds
eyyno0 = exco epyy epsno0 ;
trac mail0 eyyno0 titr 'Epsilon yy' ;
***********************************************
******** contrainte équivalentes Von-Mises dans l'élément
vmis0 = vmis mod0 sig0 ;
titr 'contrainte équivalentes Von-Mises' ;
trac mod0 vmis0 ;
******** moyenne aux noeuds
vmisno0 = chan chpo mod0 vmis0 ;
trac mail0 vmisno0 ;
trac mail0 vmisno0 (cont mail0) ;
******** contrainte équivalentes Tresca dans l'élément
tres0 = tres mod0 sig0 ;
titr 'contrainte équivalentes Tresca' ;
trac mod0 tres0 ;
******** moyenne aux noeuds
tresno0 = chan chpo mod0 tres0 ;
trac mail0 tresno0 (cont mail0) ;
***** avec les mêmes isovaleurs
isov0 = prog 0 75 150 225 300 375 450 525 600 675 750 825 900 ;
isov0 = prog 0 50 100 150 200 250 300 350 400 450 500 550 600 650
                  700 750 800 850 900 ;
trac mail0 vmisno0 (cont mail0) isov0 
 titr 'contrainte équivalentes Von-Mises';
trac mail0 tresno0 (cont mail0) isov0 
 titr 'contrainte équivalentes Tresca';
*
*
*
*********1*********2*********3*********4*********5*********6*********7
******** convergence de la solution en déplacement (en micronmètre)
*********1*********2*********3*********4*********5*********6*********7
** élément TRI6
* n              5     8    10    14    20    30    40
nnoe6 = prog   161   389   605  1169  2409  5465  9697 ;
u6 = prog     1.72  1.75  1.91  1.91  2.03  2.07  2.30 ;
v6 = prog    17.49 19.77 20.60 22.25 23.78 25.13 26.83 ;
un6 = evol bleu manu nnoe6 u6 ;
vn6 = evol bleu manu nnoe6 v6 ;
*********1*********2*********3*********4*********5*********6*********7
** élément TRI3
* n           2     3     4     5     8    10    12    13    14
*            15    40    80   100
nnoe3 = prog 13    19    30    46   106   162   232   268   307
            357  2465  9771 15236 ;
u3 = prog  1.49  0.89  1.55  1.57  1.59  1.76  1.85  1.62  1.70
           1.73  2.12  2.16  2.28 ;
v3 = prog  9.47  9.91 11.97 12.85 15.11 15.78 16.61 17.30 17.55
          17.96 22.02 24.74 25.68 ;
un3 = evol rouge manu nnoe3 u3 ;
vn3 = evol rouge manu nnoe3 v3 ;
*********1*********2*********3*********4*********5*********6*********7
** élément QUA4
* n            5     8    12    20    27    35   40     50    55    70   
nnoe4 = prog  42   103   206   523   937  1634  2101  2848  4312  7256 ;
u4 = prog   1.99  1.88  2.06  2.17  2.13  2.44  2.42  2.53  2.55  2.72 ;
v4 = prog  14.11 14.80 16.92 18.04 20.71 20.23 20.67 21.88 22.31 23.46 ;
un4 = evol vert manu nnoe4 u4 ;
vn4 = evol vert manu nnoe4 v4 ;
*********1*********2*********3*********4*********5*********6*********7
** élément QUA8
* n             4     7    10    24    27    35    40    55
nnoe8 = prog   86   230   428  2255  2790  4848  6241 13700 ;
u8 = prog    1.66  2.14  2.40  2.47  2.65  2.58  2.58  3.91 ;
v8 = prog   77.24 19.18 19.72 22.76 24.77 24.38 24.82 27.79 ;
un8 = evol jaune manu nnoe8 u8 ;
vn8 = evol jaune manu nnoe8 v8 ;
*********1*********2*********3*********4*********5*********6*********7
t = table ; tt = table ;
t.1 = 'MARQ CARR' ; tt.1 = 'TRI6' ;
t.2 = 'MARQ LOSA' ; tt.2 = 'TRI3' ;
t.3 = 'MARQ CROI' ; tt.3 = 'QUA4' ;
t.4 = 'MARQ CROI' ; tt.4 = 'QUA8' ;
t.titre = tt ;
dess (un6 et un3 et un4 et un8) lege t xbor 0.0 16000.0 
 titx 'Nb.noeuds' tity '-Ux(P4) en micron m'
 titr 'Convergence solution' ;
dess (un6 et un3 et un4 et un8) lege t xbor 0.0 16000.0 ybor 0.0 3.0 
 titx 'Nb.noeuds' tity '-Ux(P4) en micron m'
 titr 'Convergence solution' ;
dess (vn6 et vn3 et vn4 et vn8) lege t xbor 0.0 16000.0
 titx 'Nb.noeuds' tity '-Uy(P4) en micron m'
 titr 'Convergence solution' ;
dess (vn6 et vn3 et vn4 et vn8) lege t xbor 0.0 16000.0 ybor 0.0 30.0 
 titx 'Nb.noeuds' tity '-Uy(P4) en micron m'
 titr 'Convergence solution' ;
fin ;