C CQ4GR1    SOURCE    PV090527  26/04/30    21:15:23     12529          

      SUBROUTINE CQ4GR1(XE,RG,IPOIN1,IDISS,IXMATR,IPMINT,iel)
C=======================================================================
C
C      CALCULE LA MATRICE DE RAIDEUR LIEE A LA VARIATION DE PRESSION DUE
C      AU MOUVEMENT SUIVANT VECZ DANS UN CHAMP DE PESANTEUR
c      élément coq4
C
C      cette matrice est non symetrique mais combinée avec la mattrice kp
c      donne une matrice symetrique c'est pourquoi on calcule les matrices
c      symetriques et dissymetriques
c
c
C
C   ENTREE
C      XE(4,4)=COODONNEES DE L ELEMENT
C      RG     =MASSE VOLUMIQUE * ACCELERATION DE GRAVITE
C      IPOIN1 =VECTEUR(POINT) DEFINISSANT LE SENS DE GRAVITE
*      IDISS : 0 ---> MATRICE SYMETRIQUE
*              1 ---> MATRICE DISSYMETRIQUE
*     IXMATR : pointeur du segment xmatri
*      ipmint:  pointeur sur le segment d'integration
C   TRAVAIL
C      XEL(4,4)   =COORDONNEES LOCALES DE L ELEMENT
C      BPSS(3,3)  =MATRICE DE PASSAGE
C      VECN(3)    =VECTEUR DEFINISSANT LA NORMALE
C      VECZ(3)    =VECTEUR DEFINISSANT LE SENS DE GRAVITE
*      rel    : matrice masse (repeur local) avec rho=1. Elle sert à
*              calculer les termes de la matrice de rigidité qui sont
*              similaires.
C
*   ON CALCULE :
C      RE(24,24)=MATRICE DE RAIDEUR REPERE GLOBAL (dans le segment xmatri)
C
C      JUILLET 95  I. POLITOPOULOS
C
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMRIGID
-INC SMINTE
-INC SMMODEL

      DIMENSION XE(3,4),XEL(3,4),BPSS(3,3),XRE(24,24)
      DIMENSION VECZ(3),VECN(3),RHOMAT(6,6), SHPWRK(6,4)
      DIMENSION BGENE(6,24),REL(24,24)

      MINTE=IPMINT
C
C     COORDONNEES LOCALES
C
      CALL CQ4LOC(XE,XEL,BPSS,IERT,0)
C
C     MISE A 0 DE LA MATRICE
C
      XMATRI = IXMATR

C
C     CONSTRUCTION DE LA NORMALE
**** il serait plus correct d'utiliser celle calculée dans cq4loc
**** mais elle n'est pas en sortie et je suis paresseux
C
      XG1 = XE(1,2) - XE(1,1)
      YG1 = XE(2,2) - XE(2,1)
      ZG1 = XE(3,2) - XE(3,1)
      XG2 = XE(1,3) - XE(1,1)
      YG2 = XE(2,3) - XE(2,1)
      ZG2 = XE(3,3) - XE(3,1)
      VECN(1) = YG1*ZG2 - ZG1*YG2
      VECN(2) = ZG1*XG2 - XG1*ZG2
      VECN(3) = XG1*YG2 - YG1*XG2
      XNORM = VECN(1)*VECN(1) + VECN(2)*VECN(2) + VECN(3)*VECN(3)
      XNORM = SQRT(XNORM)
        DO 10 IA=1,3
          VECN(IA) = VECN(IA)/XNORM
  10    CONTINUE

C
C    DETERMINATION ET NORMALISATION DE VECZ
C
      IF (IPOIN1.NE.0) THEN
       VECZ(1) = XCOOR((IPOIN1-1)*(IDIM+1) + 1)
       VECZ(2) = XCOOR((IPOIN1-1)*(IDIM+1) + 2)
       VECZ(3) = XCOOR((IPOIN1-1)*(IDIM+1) + 3)
      ELSE
       VECZ(1) = VECN(1)
       VECZ(2) = VECN(2)
       VECZ(3) = VECN(3)
      ENDIF
      XNORM = VECZ(1)*VECZ(1) + VECZ(2)*VECZ(2) + VECZ(3)*VECZ(3)
      XNORM = SQRT(XNORM)
        DO 20 IA=1,3
          VECZ(IA) = VECZ(IA)/XNORM
  20    CONTINUE


c
c   construction de la matrice masse dans le repere local avec rho=1
c
         CALL ZERO(RHOMAT,6,6)
         RHOMAT( 1, 1)=1.D0
         RHOMAT( 2, 2)=1.D0
         RHOMAT( 3, 3)=1.D0
         CALL ZERO(REL,24,24)
         DO 21 IGAU=1,4
            CALL NCOQ4(IGAU,XEL,SHPTOT,SHPWRK,BGENE,DJAC,IERT)
C           IERT=1 JACOBIANO=<0
            IF(IERT.EQ.1) CALL ERREUR(664)
            DJAC=DJAC*POIGAU(IGAU)
            LRE = 24
            CALL BDBST(BGENE,DJAC,RHOMAT,LRE,6,REL)
 21      CONTINUE


C
C     TERMES DE LA MATRICE DANS LE REPERE GLOBAL
C
      DO 100 IA=1,19,6
       DO 110 IC= 1,3
        IAA = IA + IC - 1
         DO 120 IB=1,19,6
          DO 130 ID = 1,3
           IBB = IB + ID -1
              XRE(IBB  ,IAA) =REL(IB,IA)*VECN(ID)*VECZ(IC)*RG
  130     CONTINUE
  120   CONTINUE
  110  CONTINUE
  100 CONTINUE

*  symetrisation de la matrice
      DO 200 IA=1,19,6
       DO 210 IC= 1,3
        IAA = IA + IC - 1
         DO 220 IB=1,19,6
          DO 230 ID = 1,3
           IBB = IB + ID -1
           IF (IDISS.EQ.0) THEN
            RE(IAA,IBB,iel) = 0.5D0*(XRE(IAA, IBB) + XRE(IBB,IAA))
           ELSE
             RE(IAA,IBB,iel) = XRE(IAA,IBB)
           ENDIF
  230     CONTINUE
  220   CONTINUE
  210  CONTINUE
  200 CONTINUE


      RETURN
      END