      SUBROUTINE etcap1d(coor,car,iopt,u,alpha,theta,sigma)
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C BUT: CALCUL DES CONTRAINTES DE L ELEMENT TRIA AP1D
C ---  
c in : coor(noe,ndim) : coordonnees R(3), Z(3) des 3 sommets
c      car, iopt      : caracteristiques des materiaux
c      U(ndim,noe): deplacements U_r et U_z aux 3 noeuds (sommets)
C      alpha(3) : coef. dilatation thermique: radial, axial, tangentiel
C      theta(3)    : theta aux 3 sommets
C out: 
C      SIGMA(4) :  S_rr, S_rz, S_zz, S_tt elastiques
C .....................................................................
      implicit none
      integer NBPOLY,NPI,NPS, IP(6), iopt, ndim, i
      PARAMETER (ndim=2, NBPOLY=3, NPI=1, NPS=1)
      DOUBLE PRECISION POIDS(NPI),P13(NBPOLY,NPI),DP13(2,NBPOLY,NPI),
     +                 POIDEL(npi),F1(npi),F2(npi),DFM1(4,4),
     +                 DFM1DP(2,3,3),XYNPI(2,NPI),alpha(3),theta(3),
     +                 COOR(3,2),car(10),elas(10),U(ndim,3),r,SIGMA(4),
     +                 estrain_zz,estrain_rr, estrain_tt, estrain_rz,
     +                 tstrain_zz,tstrain_rr, tstrain_tt,alpha_z,
     +                 alpha_r,alpha_t,p1,p2,p3,dp1dr,dp2dr,dp3dr,
     +                 dp1dz,dp2dz,dp3dz
C     -- XYNPI: coordonnees pt. int. numeriques (element reference)
      DATA XYNPI /  .33333333333333333D+00,  .33333333333333333D+00 /
C     -- POIDS: poids du schema d'integration numerique.
      DATA POIDS/ .50000000000000000D+00 /
C     -- Valeurs des Polynomes de base aux pt. int. numerique.
      DATA P13 / .33333333333333333D+00,  .33333333333333333D+00,
     +           .33333333333333333D+00 /
C     -- Valeurs Derive'es des Poly. de base aux pt int. numerique.
      DATA DP13/ -1.0D+00, -1.0D+00, 1.0D+00, .0D+00, .0D+00, 1.0D+00/
C
C     CALCUL DE F1,F2,DFM1,DFM1DP
C     -----------------------
      CALL E1AP1D(NBPOLY,NPI,POIDS,P13,DP13,ip,f1,f2,DFM1dp,poidel,
     +            dfm1,COOR)
C
      call hookax(iopt,car,ELAS)
C
C     Contriantes [SIG] = [Elas] * [STRAIN_elas - STRAIN_ther] + SIG_0
C

C     1) --- Deformations Elastiques
c                                                                  [ U_r(1) ]
c  [  Strain_zz  ]    [   0      0      0    dp1dz  dp2dz  dp3dz ] [        ]
c  [             ]    [                                          ] [ U_r(2) ]
c  [  Strain_rr  ]    [ dp1dr  dp2dr  dp3dr    0      0      0   ] [        ]
c  [             ]    [                                          ] [ U_r(3) ]
c  [  Strain_tt  ] := [   p1     p2     p3                       ] [        ]
c  [             ]    [  ----   ----   ----    0      0      0   ] [ U_z(1) ]
c  [ 2*Strain_rz ]    [    r      r      r                       ] [        ]
c                     [                                          ] [ U_z(2) ]
c                     [ dp1dz  dp2dz  dp3dz  dp1dr  dp2dr  dp3dr ] [        ]
c                                                                  [ U_z(3) ]
c
      i  = 1
      r = F1(i)
      p1 = P13(1,i)
      p2 = P13(2,i)
      p3 = P13(3,i)
      dp1dr = DFM1DP(1,1,i)
      dp2dr = DFM1DP(1,2,i)
      dp3dr = DFM1DP(1,3,i)
      dp1dz = DFM1DP(2,1,i)
      dp2dz = DFM1DP(2,2,i)
      dp3dz = DFM1DP(2,3,i)

      estrain_zz = dp1dz*U(2,1) + dp2dz*U(2,2) + dp3dz*U(2,3)
      estrain_rr = dp1dr*U(1,1) + dp2dr*U(1,2) + dp3dr*U(1,3)
      estrain_tt = ( p1*U(1,1) + p2*U(1,2) + p3*U(1,3) )/r
      estrain_rz = .5*( dp1dz*U(1,1) + dp2dz*U(1,2) + dp3dz*U(1,3)
     +                + dp1dr*U(2,1) + dp2dr*U(2,2) + dp3dr*U(2,3) )
C
C     2) --- Deformations THermiques =  (ALPHA) * (P) * THETA
c         
c  [  Strain_zz ]  [alpha_z]
c  [            ]  [       ]
c  [  Strain_rr ]  [alpha_r]                  [ theta(1) ]
c  [            ]:=[       ] * [ p1 p2 p3 ] * [ theta(2) ]
c  [  Strain_tt ]  [alpha_t]                  [ theta(3) ]
c  [            ]  [       ] 
c  [ 2*Strain_rz]  [  0    ] 
c
      alpha_r = alpha(1)
      alpha_z = alpha(2)
      alpha_t = alpha(3)
      tstrain_zz = alpha_z *(p1*theta(1) + p2*theta(2) + p3*theta(3))
      tstrain_rr = alpha_r *(p1*theta(1) + p2*theta(2) + p3*theta(3))
      tstrain_tt = alpha_t *(p1*theta(1) + p2*theta(2) + p3*theta(3))

C     3) --- Contraintes
c
c  [  Sig_zz  ]    [ elas(1)  elas(2)  elas(4)   elas(7) ] [  Strain_zz  ]
c  [          ]    [                                     ] [             ]
c  [  Sig_rr  ]    [ elas(2)  elas(3)  elas(5)   elas(8) ] [  Strain_rr  ]
c  [          ] := [                                     ] [             ]
c  [  Sig_tt  ]    [ elas(4)  elas(5)  elas(6)   elas(9) ] [  Strain_tt  ]
c  [          ]    [                                     ] [             ]
c  [  Sig_rz  ]    [ elas(7)  elas(8)  elas(9)  elas(10) ] [ 2*Strain_rz ]
c
c    Rangement SIGMA(4): S_rr, S_rz, S_zz, S_tt

      sigma(3) = elas(1)*(estrain_zz - tstrain_zz)
     +         + elas(2)*(estrain_rr - tstrain_rr)
     +         + elas(4)*(estrain_tt - tstrain_tt)
     +         + 2*elas(7)*estrain_rz
      sigma(1) = elas(2)*(estrain_zz - tstrain_zz)
     +         + elas(3)*(estrain_rr - tstrain_rr)
     +         + elas(5)*(estrain_tt - tstrain_tt)
     +         + 2*elas(8)*estrain_rz
      sigma(4) = elas(4)*(estrain_zz - tstrain_zz)
     +         + elas(5)*(estrain_rr - tstrain_rr)
     +         + elas(6)*(estrain_tt - tstrain_tt)
     +         + 2*elas(9)*estrain_rz
      sigma(2) = elas(7)*(estrain_zz - tstrain_zz)
     +         + elas(8)*(estrain_rr - tstrain_rr)
     +         + elas(9)*(estrain_tt - tstrain_tt)
     +         + 2*elas(10)*estrain_rz

      return
      END
