***** dimension 3, element tétraèdre à 4 noeuds et 3 ddl par noeud
*****              les déplacements suivant chacune des 3 directions
opti dime 3 elem tet4 echo 0 ;
********** distance approximative entre 2 noeuds sur le contour
dens 1.50 ;
***************************** les points
p1 =  0.0   0.0    0.0 ;
p2 =  0.0  15.0    0.0 ;
p3 =  0.0  15.0    1.844 ;
c1 =  0.0  13.0    1.844 ;
p4 =  0.0  13.04   3.844 ;
p5 =  0.0   6.061  3.983 ;
c2 =  0.0   6.141  7.982 ;
p6 =  0.0   2.142  7.902 ;
p7 =  0.0   1.718 29.020 ;
c3 =  0.0   0.718 29.0 ;
p8 =  0.0   0.718 30.0 ;
p9 =  0.0   0.0   30.0 ;
***************************** les lignes
l12 = droi p1 p2 ;
l23 = droi p2 p3 ;
l34 = cerc p3 c1 p4 ;
l45 = droi p4 p5 ;
l56 = cerc p5 c2 p6 ;
l67 = droi p6 p7 ;
l78 = cerc p7 c3 p8 ;
l89 = droi p8 p9 ;
************************* le demi contour
cont1= l12 et l23 et l34 et l45 et l56 et l67 et l78 et l89 ;
************************* on symétrise
cont2 = cont1 syme droi p1 p9 ;
************************* assemble les 2 parties du contour
cont0 = cont1 et cont2 ;
************************* élimine les noeuds doubles
elim cont0 ; trac cont0 ;
************************* fabrique le maillage surfacique
sur0 = surf cont0 plan ; trac sur0 ;
************************* on fabrique le maillage volumique
************************* par translation de la surface
************************* suivant le vecteur vec1
lo0 = 50.0 ;
vec1 = lo0 0.0 0.0 ;
mail0 = sur0 volu tran vec1 ; 
trac mail0 ; trac cach mail0 ;
nbno0 = nbno mail0 ; mess 'Nombre de noeuds :' nbno0 ;
*
********************* 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
cl0 = bloq depl sur0 ;
*
******** recherche du noeud (pf0) le plus proche des coordonnées fournies
pf0 = mail0 poin proc (lo0 0.0 0.0) ;
*pf0 = mail0 poin proc (lo0 5.0 0.0) ;
******** on applique une force ponctuelle fo0 au noeud pf0 
fo0 = 1000.0 ;
vfo0 = 0.0 0.0 fo0 ;
for0 = forc vfo0 pf0 ;
******** pour le dessin
amp0 = -4.0e-2 ;
vect0 = vect for0 amp0 fx fy fz vert ;
trac cach vect0 mail0 titr 'Chargements' ;
*
********************* 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' ;
trac cach (def0 et def1) titre 'Déformée amplifiée' ;
def1 = defo res0 mail0 1.0 rouge ;
trac cach (def0 et def1) titre 'Déformée à l échelle de la structure' ;
******** 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 le dessin on met la même amplitude
***** que pour les forces extérieurs imposées
vreac0 = vect reac0 amp0 fx fy fz vert ;
trac vreac0 mail0 titr 'Action aux liaisons';
trac cach 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 ;
*
*********************** 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 ;
***** on exTRAIT la coMPOSANTE de nom sm.. 
***** des tenseurs de contrainte (1 tenseur par élément)
***** contenus dans sig0 pour la tracer
sxx0 = exco smxx sig0 ;
trac mod0 sxx0 titr 'Sigma xx' ;
sxy0 = exco smxy sig0 ;
trac mod0 sxy0 titr 'Sigma xy' ;
syy0 = exco smyy sig0 ;
trac mod0 syy0 titr 'Sigma yy' ;
sxz0 = exco smxz sig0 ;
trac mod0 sxz0 titr 'Sigma xz' ;
syz0 = exco smyz sig0 ;
trac mod0 syz0 titr 'Sigma yz' ;
szz0 = exco smzz sig0 ;
trac mod0 szz0 titr 'Sigma zz' ;
****** calcul des contraintes aux noeuds
****** à partir des contraintes dans les éléments
****** signo0 contiendra les valeurs du tenseur des contraintes
****** pour chacun des noeuds du maillage
signo0 = chan chpo mod0 sig0 ;
******** contrainte équivalentes Von-Mises dans les éléments
vmis0 = vmis mod0 sig0 ;
titr 'contrainte équivalentes Von-Mises' ;
trac mod0 vmis0 ;
******** on coupe la structure suivant un plan pour y voir
******** les contrainte équivalentes Von-Mises
po3 = lo0 0.0 0.0 ;
trac mod0 vmis0 coup p1 p9 po3 ;
******** on regarde les actions aux liaisons vu de oeil1
oeil1 = 0.0 ((-1.0)*lo0) 0.0 ;
trac cach oeil1 vreac0 mail0 titr 'Action aux liaisons' ;
fin ;