opti dime 2 elem tri3 mode plan cont echo 1 dens 1.0 ;
******** 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
*
********************* géométrie
***** les points définissant le contour extérieur
p1 =  0.0  0.0 ;
p2 =  4.0  0.0 ;
p3 =  8.0  0.0 ;
c1 =  6.0 -2.0 ;
p4 = 12.0  0.0 ;
p5 = 10.0  5.0 ;
p6 =  0.0  5.0 ;
***** les lignes définissant le contour extérieur
l1 = droi p1 p2 ;
l2 = cerc p2 c1 p3 ;
l3 = droi p3 p4 ;
l4 = droi p4 p5 ;
l5 = droi p5 p6 ;
l6 = droi p6 p1 ;
***** le contour extérieur
cont0 = l1 et l2 et l3 et l4 et l5 et l6 ;
trac cont0 ;
******** maillage
mail0 = surf cont0 plan ;
trac mail0 titr 'Maillage' ;
nbno0 = nbno mail0 ; mess 'Nombre de noeuds :' nbno0 ;
nbel0 = nbel mail0 ; mess 'Nombre d éléments :' nbel0 ;
*
********************* modèle (homogène et linéaire sous entendu)
mod0 = mode mail0 mecanique elastique isotrope ;
mat0 = mate mod0 youn 210000.0 nu 0.30 ;
*
********************* conditions limites
******** blocage du déplacement suivant x et y sur l1
******** blocage du déplacement suivant y sur l3
cl1 = bloq depla l1 ;
cl2 = bloq uy l3 ;
cl0 = cl1 et cl2 ;
*
********************* chargement ponctuel
********** force ponctuelle de 500.0 suivant x en p5
******** pour le calcul 
for1 = forc fx 500.0 p5 ;
******** pour le dessin (amp0 est l'amplitude de représentation)
amp0 = 1.0e-2 ;
vect1 = vect for1 amp0 fx fy vert ;
******** force répartie constante de 100.0 suivant y sur l5
vfr2 = 0.0 -100.0 ;
******** pour le calcul
for2 = fsur mass mod0 l5 vfr2 ;
******** pour le dessin (avec la même amplitude)
vect2 = vect for2 amp0 fx fy vert ;
trac (vect2 et vect1) mail0 titr 'Chargements' ;
******** le second membre total pour le calcul
for0 = for1 et for2 ;
*
********************* matrice de raideur
rig0 = rigi mod0 mat0 ;
******** prise en compte des conditions limites dans la matrice de raideur
rig0 = rig0 et cl0 ;
*
********************* résolution de [K][q]=[F] avec
******** [K] = rig0 ; [F] = for0 et [q] = res0 contient les ddl (u et v)
******** en chacun des noeuds du maillage
res0 = reso rig0 for0 ;
*
********************* déformée
def0 = defo res0 mail0 0.0 ;
def1 = defo res0 mail0 rouge ;
trac (def0 et def1) titre 'Déformée amplifiée' ;
def1 = defo res0 mail0 1.0 rouge ;
trac (def0 et def1) titre 'Déformée à l échelle de la structure' ;
*
******** valeur déplacement en P5
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 '' ;
*
******** valeur déplacement au poinT le plus procHE de (5,2.5)
mess 'x,y désiré : ' 5.0 2.5 ;
p7 = mail0 poin proc (5.0 2.50) ;
x7 = coor 1 p7 ; y7 = coor 2 p7 ;
mess 'x,y du point P7 le plus proche : 'x7 y7 ;
ux7 = extr res0 ux p7 ;
mess 'Déplacement suivant x en P7 : 'ux7 ;
uy7 = extr res0 uy p7 ;
mess 'Déplacement suivant y en P7 : 'uy7 ; mess '' ;
*
******** déplacement aux noeuds de l4
****** on créé une courbe (avec evolUTION) 
****** à partir du chAMP PAR poINT contenu dans res0
****** et ayant le nom ux et uniquement pour les points de l4
uxl4 = evol bleu chpo res0 ux l4 ;
****** idem pour uy
uyl4 = evol rouge chpo res0 uy l4 ;
****** pour la légende
t = table ; tt = table ;
t.1 = 'MARQ CARR' ; tt.1 = 'Ux' ;
t.2 = 'MARQ LOSA' ; tt.2 = 'Uy' ;
t.titre = tt ;
****** on dessINE les 2 courbes
dess (uxl4 et uyl4) lege t  
 titx 'ligne l4' tity 'Ux et Uy'
 titr 'Ux et Uy sur l4' ;
*
******** action aux liaisons
***** reac0 contient les forces nodales pour les noeuds où l'on a imposé
***** les conditions aux limites
reac0 = reac rig0 res0 ;
***** pour visualiser le contenu d'une variable créé
*list reac0 ;
***** pour le dessin on met la même amplitude
***** que pour les forces extérieurs imposées
vreac0 = vect reac0 amp0 fx fy vert ;
trac vreac0 mail0 titr 'Action aux liaisons';
***** on veut avoir la force globale (la resuLTANTE de ces forces nodales)
fg0 = resu reac0 ;
mess 'Force globale :' ; list fg0 ;
***** on veut avoir la force globale sur chacun des contours l1 et l3
***** (reduIRE)
reac1 = redu reac0 l1 ; fg1 = resu reac1 ;
mess 'Force sur l1 :' ; list fg1 ;
reac1 = redu reac0 l3 ; fg1 = resu reac1 ;
mess 'Force sur l3 :' ; list fg1 ;
***** (les anciennes variables reac1, fg1 sont écrasées)
*
*********************** calcul des contraintes dans l'élément
****** sig0 contiendra les valeurs du tenseur des contraintes
****** pour chacun des éléments du maillage
sig0 = sigm mod0 mat0 res0 ;
******** tracé de sIGmAxx sur les elements
***** on exTRAIT la coMPOSANTE de nom sIGmAXX des tenseurs contenus dans sig0
sxx0 = exco smxx sig0 ;
trac mod0 sxx0 titr 'Sigma xx' ;
******** tracé de sIGmAxy sur les elements
sxy0 = exco smxy sig0 ;
trac mod0 sxy0 titr 'Sigma xy' ;
******** tracé de sIGmAyy sur les elements
syy0 = exco smyy sig0 ;
trac mod0 syy0 titr 'Sigma yy' ;
******** pour faire beau on moyenne les contraintes aux noeuds
****** signo0 contiendra les valeurs du tenseur des contraintes
****** pour chacun des noeuds du maillage
signo0 = chan chpo mod0 sig0 ;
******** tracé de sIGmAxx moyenné aux noeuds
sxxno0 = exco smxx signo0 ;
trac mail0 sxxno0 titr 'Sigma xx' ;
******** tracé de sigmaxy moyenné aux noeuds
sxyno0 = exco smxy signo0 ;
trac mail0 sxyno0 titr 'Sigma xy' ;
******** trace de sigmayy moyenné aux noeuds
syyno0 = exco smyy signo0 ;
trac mail0 syyno0 titr 'Sigma yy' ;
******** rem : mail0 contient des renseignements aux noeuds
******** alors que mod0 contient des renseignements dans les éléments
* sur les graphes précédents, une couleur déterminée
* ne correspond pas toujours à la même valeur d'où le tracé suivant
* list0 ne doit pas comporter plus de 14 valeurs !
* de plus le tracé ne se fait pas sur le maillage initial mais sur la déformée (non amplifiée)
list0 = prog -250.0 -200.0 -150.0 -100.0  -50.0     0.0   50.0 
                    100.0   150.0  200.0   250.0 300.0 350.0 400.0 ;
trac list0 sxxno0 def1 titr 'Sigma xx' ;
trac list0 sxyno0 def1 titr 'Sigma xy' ;
trac list0 syyno0 def1 titr 'Sigma yy' ;
* rem : l'ordre des opérateurs n'est pas important
******** représentation des contraintes sur l2
sxxl2 = evol rouge chpo sxxno0 l2 ;
sxyl2 = evol bleu chpo sxyno0 l2 ;
syyl2 = evol vert chpo syyno0 l2 ;
t = table ; tt = table ;
t.1 = 'MARQ CARR' ; tt.1 = 's xx' ;
t.2 = 'MARQ LOSA' ; tt.2 = 's xy' ;
t.3 = 'MARQ CROI' ; tt.3 = 's yy' ;
t.titre = tt ;
dess (sxxl2 et sxyl2 et syyl2) lege t  
 titx 'ligne l2' tity 'sigma'
 titr 'sxx, sxy et syy sur l2' ;
*
******** 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' ;
*
******** densité d'énergie de déformation dans l'élément
****** ener0 est un scalaire calculé sur chaque élément
ener0 = ener mod0 sig0 eps0 ;
titr 'Densité surfacique d énergie (J/mm2)' ;
trac mod0 ener0 ;
******** moyenne des densités d'énergie aux noeuds
enerno0 = chan chpo mod0 ener0 ;
trac mail0 enerno0 ;
*
******** 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 cont0 ;
*
******** 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 cont0 ;
******** comparaison de Tresca et Von Mises 
* avec les mêmes isovaleurs
list0 = prog 0.0 50.0 100.0 150.0 200.0 250.0 300.0
                       350.0 400.0 450.0 500.0 550.0 600.0 650.0 ;
trac list0 mail0 vmisno0 cont0
       titr 'contrainte équivalentes Von-Mises' ;
trac list0 mail0 tresno0 cont0
       titr 'contrainte équivalentes Tresca' ;
trac mail0 (2.0*(abs(tresno0-vmisno0))/(tresno0+vmisno0)) cont0 
       titr 'erreur relative (en %) entre Tresca et Von-Mises';
*
fin ;