      SUBROUTINE ETRAP1D(coor,car,iopt,AE)
C  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C  BUT  : CALCUL DE LA MATRICE DE RAIDEUR DE L ELEMENT TRIA AP1D
C  ---  
c in : coor(noe,ndim) : coordonnees R(3), Z(3) des 3 sommets.
c      car            : caracteristiques des materiaux
c           if(iopt .lt. 3) 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: ae            : matrice triangulaire sup
c ...............................................................
c  programmeur : modulef
c ...............................................................
      implicit none
      integer ndim, NBPOLY, NPI, IP(6), l, iopt, i
      PARAMETER (ndim=2 , NBPOLY=3 , NPI=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),ae(21),car(10),
     +                 COOR(3,2),elas(10),AUX,A11(1),A12(2),A13(2),
     +                 A22(2,2),A23(2,2),A33(3),G1(3,3),G2(3),G3(3,2),
     +                 G5(2,3),G4(3,2)
C     -- XYNPI: coordonnees pt. int. numeriques (element reference)
      DATA XYNPI /  .33333333333333333D+00,  .33333333333333333D+00 /
C     -- POIDSG: 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
      CALL e1ap1d(nbpoly,npi,poids,p13,dp13,ip,f1,f2,dfm1dp,poidel,
     +            dfm1,coor)
C
      call hookax(iopt,car,elas)
C
C     -- INITIALISATION DE AE
      do 1 i = 1 , 21
        ae(i) = 0.d0
 1    continue
C
C     CONSTRUCTION DES MATRICES   (A11 , A12 , A13)
C     -------------------------   (      A22 , A23)  =  (TD * E * D)
C                                 (            A33)
      DO 2 L=1,NPI
         AUX      = POIDEL(L) / F1(L)
         A11(1)   = AUX * ELAS(6) / F1(L)
         A12(1)   = AUX * ELAS(5)
         A12(2)   = AUX * ELAS(9)
         A13(1)   = AUX * ELAS(9)
         A13(2)   = AUX * ELAS(4)
         AUX      = POIDEL(L)
         A22(1,1) = AUX * ELAS(3)
         A22(2,1) = AUX*ELAS(8)
         A22(1,2) = A22(2,1)
         A22(2,2) = AUX * ELAS(10)
         A23(1,1) = A22(2,1)
         A23(2,1) = A22(2,2)
         A23(1,2) = AUX * ELAS(2)
         A23(2,2) = AUX * ELAS(7)
         A33(1)   = A22(2,2)
         A33(2)   = A23(2,2)
         A33(3)   = AUX * ELAS(1)
C
C        CALCUL DE G1
C        ------------
         CALL AB0D(1,1,NBPOLY,A11,P13(1,L),G2)
         CALL AB1D(1,2,NBPOLY,A12,DFM1DP(1,1,L),G2)
C
         CALL TAB0D(NBPOLY,1,NBPOLY,P13(1,L),G2,G1)
C
         CALL TAB0D(2,1,NBPOLY,A12,P13(1,L),G5)
         CALL AB1D(2,2,NBPOLY,A22,DFM1DP(1,1,L),G5)
C
         CALL TAB1D(NBPOLY,2,NBPOLY,DFM1DP(1,1,L),G5,G1)
C        ON PLONGE G1 DANS AE
         CALL PLONAD(-1,NBPOLY,NBPOLY,1,1,IP,G1,G1,AE)
C
C        CALCUL DE G2
C        ------------
         CALL TAB0D(NBPOLY,1,2,P13(1,L),A13,G3)
         CALL TAB1D(NBPOLY,2,2,DFM1DP(1,1,L),A23,G3)
C
         CALL AB0D(NBPOLY,2,NBPOLY,G3,DFM1DP(1,1,L),G1)
C        ON PLONGE G2 DANS AE
         CALL PLONAD(-1,NBPOLY,NBPOLY,1,2,IP,G1,G1,AE)
C
C        CALCUL DE G3
C        ------------
         CALL TABAXD(NBPOLY,2,DFM1DP(1,1,L),A33,G1,G4,1)
C        ON PLONGE G3 DANS AE
         CALL PLONAD(1,NBPOLY,NBPOLY,2,2,IP,G1,G1,AE)
 2    CONTINUE
      END
