champgeoquad "ChampGeo" "cy4t-1000.geoquad"
scalaire "AlphaMin" .01002
scalaire "AlphaMax" .01014
scalaire "MaxIt" 11
scalaire "Epais" 0.001
scalaire "p0" 0.625-0.5625*x*x+0.0625*x*x*x*x
scalaire "p1" f(Epais)=(-0.0543421+0.06521186*x*x)*sqrt(Epais)
scalaire "p2" f(Epais)=(-0.51078287 +0.18427874*x*x)*Epais
scalaire "p3" f(Epais)=(-1.8902477-0.11709027*x*x)*Epais*sqrt(Epais)
scalaire "p4" f(Epais)=(-7.04989265+7.88572271*x*x-1.38020833*x*x*x*x+0.04921875*x*x*x*x*x*x-0.00050223*x*x*x*x*x*x*x*x)*Epais*Epais
scalaire "p5" f(Epais)=(4.64512637-0.54520440*x*x+0.10393140*x*x*x*x-0.00244544*x*x*x*x*x*x)*Epais*Epais*sqrt(Epais)
scalaire "a"  f(Epais)=1.27788621/sqrt(Epais) + 1.56508458*sqrt(Epais) + 0.98503729*Epais*sqrt(Epais) - 2.91823062*Epais*Epais*sqrt(Epais)
scalaire "b"  f(Epais)=1.27788621/sqrt(Epais) + 1.56508458*sqrt(Epais) + 0.98503729*Epais*sqrt(Epais) + 2.91823062*Epais*Epais*sqrt(Epais)
scalaire "C"  f(Epais)=-0.125-0.0108684*sqrt(Epais) + 0.32650413*Epais + 2.00733804*Epais*sqrt(Epais) + 0.49566176*Epais*Epais - 4.20140792*Epais*Epais*sqrt(Epais)
scalaire "S"  f(Epais)=0.125+0.69559315*sqrt(Epais) - 0.12237999*Epais - 1.43051440*Epais*sqrt(Epais) - 0.86217411*Epais*Epais - 8.78638051*Epais*Epais*sqrt(Epais)
scalaire "couchelim" f(a,b,C,S)=exp(a*(x-1))*(C*cos(b*(x-1))+S*sin(b*(x-1))) + exp(-a*(x+1))*(C*cos(b*(x+1))-S*sin(b*(x+1)))

scalaire "VTheo" f(p0,p1,p2,p3,p4,p5,couchelim)=(5.33333333333333333e-5)*(p0 + p1 + p2 + p3 + p4 + p5 + couchelim)
//
scalaire "Normal" f(Epais)=-Epais*(z*z-y*y)/((y*y+z*z)*sqrt(y*y+z*z))
scalaire "p_0y" f(Normal)=Normal*y
scalaire "p_0z" f(Normal)=Normal*z
scalaire "zero1d" 0.
vectoriel3d "zero3d" [0.,0.,0.]
empilementv3d "p_0" [zero1d,p_0y,p_0z]
//
scalaire "E1" 2.e5
scalaire "Nu1" 0.333333333
//
// le champs de deplacements initialise a zero:
//
v3dcrourav "U" [0,0,0]
//
// pour le non-conforme
//
v3dcrourav "UNc" [0,0,0]
//
//pour les reactions
//
//v3dcrourav "R" [0,0,0]
//
// les forces volumiques (si <> 0)
//vectoriel3d "F" [0,0,-8000.] 
//
// le tenseur d'elasticite (isotrope + pas de changement de repere + constant)
//
//[E1,E2,E3,Nu1,Nu2,Nu3,E1/(2.0*(1+Nu1),E2/(2.0*(1+Nu2)),E3/(2.0*(1+Nu3))] 
scalaire "G1" f(Nu1,E1)=E1/(2.0*(1+Nu1))
//
to4elasticite "E" [E1,E1,E1,Nu1,Nu1,Nu1,G1,G1,G1]
//