      SUBROUTINE etsap1d(coor,fomega,fgamma,pressi,norefs,alpha,
     +                   theta,car,iopt,be)
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C BUT: CALCUL DE SECONDS MEMBRES DE L ELEMENT TRIA AP1D
C ---  
c in : 
c   coor(noe,ndim) : coordonnees R(3), Z(3) des 3 sommets.
c   fomega(ndim,noe)        : fx, fy aux noeuds => fomega(ndim,npi)
c                               => fface(ndim=2 , npi=1)
c   fgamma(ndim,nnof*nbarete): fx, fy aux noeuds de chaque arete.
c   pressi(nnof*nbarete)     : pression aux noeuds de chaque arete
c                               => farete(ndim=2  ,npia*nbarete=6)
c   norefs(nbarete,2): norefs(i,1) = 0 si fgamma   = 0 sur arete_i
c                      norefs(i,2) = 0 si pression = 0 sur arete_i
C   theta(3)    : theta aux 3 sommets
C   alpha    : coef. dilatation thermique: radial, axial, tangentiel
c   car        : caracteristiques des materiaux
c        if(iopt .eq. 1) then
c          car(1) = young
c          car(2) = poisson
c        else
c          car(1) = E_1  (Young radial     -> E_r     )
c          car(2) = nu_1 (poisson          -> Nu_r    )
c          car(3) = E_2  (Young axial      -> E_z     )
c          car(4) = nu_2 (poisson          -> Nu_z    )
c          car(5) = E_3  (Young Tangentiel -> E_theta )
c        end if
C out: BE(6): second membre.
C  ..................................................................
      implicit none
      integer NBPOLY,NPI,IP(6),npia,i,j,i1,l,k,k1,l1,iopt,norefs(3,2),
     +        nbarete, nnof, il
      PARAMETER (NPI=1, NBPOLY=3, nbarete=3, nnof=2)
      DOUBLE PRECISION POIDS(NPI),P13(NBPOLY,NPI),DP13(2,NBPOLY,NPI),
     +                 POIDEL(npi),F1(npi),F2(npi),DFM1(4,4),
     +                 DFM1DP(2,3,3),fomega(2,3),fface(2,npi),
     +               fgamma(2,6),pressi(6),farete(2,6),xnu,ynu,alpha(3),
     +                 alpha_r,alpha_z,alpha_t,theta(3),BE(6),COOR(3,2),
     +                 car(10),D1,RR,D2PI,POIDSA(2),XA(2),D,
     +                 elas(10),A(4),G1(6),G2(6),G4(6,3),XYNPI(2,NPI),
     +                 fx,fy,xji,yji,arelon(nbarete),
     +                 p13a(2,nnof)
      DATA D2PI/ 6.2831853071795862D0 /
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     -- integration / aretes
      DATA NPIA,POIDSA,XA/2, 0.5D0,0.5D0,
     +                       0.211324865404D0,0.788675134593D0/

c   -- Valeurs aux pt d'int. num. a partir valeurs aux noeuds
c      Efforts volumiques  fomega(ndim,noe) -> fface(ndim,npi)
      do 2 i = 1 , npi
        fface(1,i) = 0.d0
        fface(2,i) = 0.d0
        do 1 j = 1 , nbpoly
          fface(1,i) =  fface(1,i) + P13(j,i) * fomega(1,j)
          fface(2,i) =  fface(2,i) + P13(j,i) * fomega(2,j)
 1      continue
 2    continue
c
c   -- Valeurs aux pt d'int. num. a partir valeurs aux noeuds
c      fgamma(ndim,(nnof*nbarete) -> farete(ndim, npia*nbarete)
c      pressi(nnof*nbarete)       -> idem
c
      do 5 k = 1 , nbarete
C       -- LONGEUR DE L ARETE, COSINUS DIRECTEURS DE LA NORMALE ext.
        J      = MOD(k,3) + 1
        XJI    = coor(J,1) - coor(k,1)
        YJI    = coor(J,2) - coor(k,2)
        ARELON(k) = DSQRT(XJI**2 + YJI**2)
        XNU    =  YJI/ARELON(k)
        YNU    = -XJI/ARELON(k)
        il = npia*(k-1)
        do 4 i = 1 , npia
          fx = 0.d0
          fy = 0.d0
          p13a(1,i) = 1.d0 - xa(i)
          p13a(2,i) = xa(i)
          do 3 j = 1 , nnof
            fx = fx + p13a(j,i) * (fgamma(1,j+il) - pressi(j+il) * XNU)
            fy = fy + p13a(j,i) * (fgamma(2,j+il) - pressi(j+il) * YNU)
 3        continue
          farete(1,il+i) = fx
          farete(2,il+i) = fy
 4      continue
 5    continue
C
c     -- CALCUL DE F1,F2,FFM1,DFM1DP
c
      CALL e1ap1d(nbpoly,npi,poids,p13,dp13,ip,f1,f2,dfm1dp,poidel,
     +            dfm1,coor)
C
      do 6 i = 1 , 6
          BE(i) = 0.D0
 6     continue
C
c     -- CONTRIBUTION DES EFFORTS SURFACIQUES 
C             FOMEGA(2,npi=3)
      DO 9 I=1,NBPOLY
         DO 8 J=1,2
            I1 = IP(I + NBPOLY * (J-1))
            D1 = 0.D0
            DO 7 L = 1,NPI
               D1 = D1 + POIDEL(L) * P13(I,L) * fface(J,L)
 7          CONTINUE
            BE(I1) = BE(I1) + D1
 8       CONTINUE
 9    CONTINUE
C
c     -- CONTRIBUTIONS DES EFFORTS SUR LES ARETES
C
      do 12 k = 1 , nbarete
        k1    = MOD(k,3) + 1
        D1    = D2PI*arelon(k)
        DO 11 L=1,NPIA
          RR = COOR(K,1) * (1-XA(L)) + COOR(K1,1) * XA(L)
          L1 = L+NPIA*(K-1)
          DO 10 J=1,2
            I=IP(K+(J-1)*NBPOLY)
            I1=IP(K1+(J-1)*NBPOLY)
            D = RR * D1 * POIDSA(L) * farete(J,L1)
            BE(I)  = BE(I)  + D * (1-XA(L))
            BE(I1) = BE(I1) + D * XA(L)
 10       CONTINUE
 11     CONTINUE
 12   CONTINUE
C
c     -- CONTRIBUTIONS DES CONTRAINTES THERMIQUES
C
      call hookax(iopt,car,elas)
c
c  [  Sig_zz  ]    [      ] [alpha_z]
c  [          ]    [      ] [       ]
c  [  Sig_rr  ]    [      ] [alpha_r]                [ theta(1) ]
c  [          ] := [ elas ]*[       ]*[ p1 p2 p3 ] * [ theta(2) ]
c  [  Sig_tt  ]    [      ] [alpha_t]           1*3  [ theta(3) ]
c  [          ]    [      ] [       ]                          3*1
c  [  Sig_rz  ]    [      ] [  0    ]
c            4*1         4*4       4*1      
C
C -- [ A ]  = [ELAS] * [ALPHA]
C       4*1       4*4       4*1
      alpha_r = alpha(1)
      alpha_z = alpha(2)
      alpha_t = alpha(3)
      A(1) = ALPHA_z*ELAS(4) + ALPHA_r*ELAS(5) + ALPHA_t*ELAS(6)
      A(2) = ALPHA_z*ELAS(2) + ALPHA_r*ELAS(3) + ALPHA_t*ELAS(5)
      A(3) = ALPHA_z*ELAS(7) + ALPHA_r*ELAS(8) + ALPHA_t*ELAS(9)
      A(4) = ALPHA_z*ELAS(1) + ALPHA_r*ELAS(2) + ALPHA_t*ELAS(4)
C
cc      CALL AZEROD(3,6,1,G4)
      do 19 I = 1 , 6
        do 19 J = 1 , 3
          G4(I,J) = 0.d0
 19   continue
C
      DO 16 L=1,NPI
        D    = POIDEL(L)
        DO 15 I=2,4
          G1(I) = D * A(I)
 15     CONTINUE
        D = D/F1(L)
        G1(1) = D * A(1)
C
C       -- TP * G1(1) + TDFM1DP * (G1(2),G1(3))
C
        CALL TAB0D(NBPOLY,1,1,P13(1,L),G1(1),G2)
        CALL TAB1D(NBPOLY,2,1,DFM1DP(1,1,L),G1(2),G2)
C
C       -- TDFM1DP * (G1(2),G1(4))
C
        CALL TAB0D(NBPOLY,2,1,DFM1DP(1,1,L),G1(3),G2(4))
C       -- G2 * P
        CALL AB1D(6,1,NBPOLY,G2,P13(1,L),G4)
 16   CONTINUE
C
C     THETA * G4  => BE = BE + G4(IP)
C
      DO 18 I=1,6
         L = IP(I)
         D = 0.D0
         DO 17 J=1,NBPOLY
            D = D + theta(J) * G4(I,J)
 17      CONTINUE
         BE(L) = BE(L) + D
 18   CONTINUE
      END
