C COQ3M1    SOURCE    CHAT      05/01/12    22:21:24     5004
      SUBROUTINE COQ3M1(XE,RHO,XEL,BPSS,BB,REWO,ILUMP)
c=======================================================================
c
c      calcule la matrice de masse de l element coq3
c
c   entree
c      xe(3,3)=coodonnees de l element
c      rho    =masse volumique * epaisseur
c      ilump = 1 si l'opérateur LUMP est appelé , 0 sinon
c   travail
c      xel(3,3)   =coordonnees locales de l element
c      bpss(3,3)  =matrice de passage
c      bb(9)      =stocke les fonctions de formes de flexion
c   sorties
c      rewo(18,18)=matrice de masse    repere local puis global
c
c      code fevrier 85   ebersolt
c
c=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XE(3,*),XEL(3,*),BPSS(3,*),REWO(18,*),BB(*)
      DIMENSION REWOLP(18,18)
      DIMENSION XX(3),YY(3)
      DATA XX/0.5D0,0.5D0,0.0D0/
      DATA YY/0.0D0,0.5D0,0.5D0/
      DATA UNSIX,UN12/.166666666666666666D0,.833333333333333333D-1/
      DATA UNTIER/.33333333333333333D0/
      DATA IZERO/0/
c
c     matrice de passage
c
      CALL VPAST(XE,BPSS)
c
c     coordonnees locales
c
      CALL VCORLC(XE,XEL,BPSS)
c
c     mise a 0 de la matrice de masse
c
      CALL ZERO(REWO,18,18)
      X21=XEL(1,2)-XEL(1,1)
      Y31=XEL(2,3)-XEL(2,1)
      SURF=X21*Y31*.5D0*RHO
c
c     termes diagonaux  de la membrane
c
      DO 100 IA=1,13,6
          REWO(IA  ,IA  )=UNSIX
          REWO(IA+1,IA+1)=UNSIX
          DO 110 IB=1,13,6
              IF(IB.EQ.IA) GOTO 110
              REWO(IA  ,IB  )=UN12
              REWO(IA+1,IB+1)=UN12
  110     CONTINUE
  100 CONTINUE
c
c     partie flexion   integration numerique
c
c     ia numero du point de gauss
c     ib numero d un noeud
c     ic numero d un noeud
c
      DO 200 IA=1,3
          CALL MFDKT(XX(IA),YY(IA),XEL,BB)
          DO 210 IB=1,3
              IBDERI=(IB-1)*6+2
              IBDEBB=(IB-1)*3
              DO 210 IC=1,3
                  ICDERI=(IC-1)*6+2
                  ICDEBB=(IC-1)*3
                  DO 220 ID=1,3
                   DO 220 IE=1,3
                    REWO(IBDERI+ID,ICDERI+IE)=REWO(IBDERI+ID,ICDERI+IE)
     1               +BB(IBDEBB+ID)*BB(ICDEBB+IE)*UNTIER
  220             CONTINUE
  210     CONTINUE
  200 CONTINUE
c
c     multiplication par rho*surf
c
      DO 300 IA=1,18
          DO 300 IB=1,18
              REWO(IA,IB)=REWO(IA,IB)*SURF
  300 CONTINUE
c
c    diagonalisation dans le cas de l'opérateur LUMP
c
c      REWO est rangé dans l'ordre i noeud x(UX UY UZ RX RY RZ) ....
c
      IF ( ILUMP .EQ. 1 ) THEN
        CALL LUMP3(REWO)
      ENDIF
c
c
c     changement de repere
c
      CALL TRANSK(REWO,BPSS,18,3,IZERO)
c
      RETURN
      END