C BGRCQ4    SOURCE    CHAT      06/03/29    21:15:59     5360
      SUBROUTINE BGRCQ4(IGAU,XEL,SHPTOT,SHP,BGR,DJAC,EXCEN,NOPLAN,IARR)
C |====================================================================|
C |                                                                    |
C |  CALCULER LA MATRICE DE GRADIENT BGR POUR COQ4 excentre ou non     |
C |     CODE SUO X.Z.                                                  |
C |====================================================================|
C |INPUT                                                               |
C |   IGAU=NUMERO DU POINT DE GAUSS                                    |
C |   NBNO=NOMBRE DE NOEUDS                                            |
C |   LRE =NOMBRE DE COLONNES DE LA MATRICE BGR                        |
C |   XEL =COORDONNEES LOCALES DE L ELEMENT                            |
C |   SHPTOT(6,NBNO,NBGAU)=FONCTIONS DE FORMES ET DERIVEES             |
C |   EXCEN=Excentrement de coque.                                     |
C |   NOPLAN=1 si element non plan,                                    |
C |         =0 sinon.                                                  |
C |   SHP=Fonctions de forme et derivees dans l'espace geometrique.    |
C |ZONE DE TRAVAIL                                                     |
C |   SHP(6,NBNO)=TABLEAU DE TRAVAIL                                   |
C |OUTPUT                                                              |
C |   DJAC=JACOBIEN                                                    |
C |   BGR(9,LRE)=MATRICE BGR                                           |
C |====================================================================|
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (LRE=24,NBNO=4)
      DIMENSION XEL(3,*),BGR(9,*),SHP(6,*),SHPTOT(6,NBNO,*)
      DIMENSION SHP7(4:5,NBNO),SHP8(4:5,NBNO)
      DATA ZER/0.D0/
      IARR=0
*
*     Matrice jacobienne:
*
      DO 110 I=1,NBNO
         SHP(1,I)=SHPTOT(1,I,IGAU)
         SHP(2,I)=SHPTOT(2,I,IGAU)
         SHP(3,I)=SHPTOT(3,I,IGAU)
 110  CONTINUE
      DXDQSI = ZER
      DXDETA = ZER
      DYDQSI = ZER
      DYDETA = ZER
      DO 120 I=1,NBNO
         DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
         DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
         DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
         DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
 120  CONTINUE
      DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
      IF(DJAC.LE. ZER) THEN
         IARR=1
         RETURN
      END IF
*
* Fonctions de forme dans l'espace geometrique
*
      DQSIDX= DYDETA/DJAC
      DQSIDY=-DXDETA/DJAC
      DETADX=-DYDQSI/DJAC
      DETADY= DXDQSI/DJAC
      DO 130 I=1,NBNO
         TT=SHP(2,I)*DQSIDX+SHP(3,I)*DETADX
         SHP(3,I)=SHP(2,I)*DQSIDY+SHP(3,I)*DETADY
         SHP(2,I)=TT
 130  CONTINUE
C
C CALCUL DE LA MATRICE BGR
C
      CALL ZERO(BGR,9,LRE)
      N1 = 1
      N2 = 2
      DX1 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2)
      DY1 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2)
      N1 = 3
      N2 = 4
      DX2 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2)
      DY2 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2)
      N1 = 1
      N2 = 4
      DX3 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2)
      DY3 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2)
      N1 = 2
      N2 = 3
      DX4 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2)
      DY4 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2)
*
      D2JAC = 2.D0 * DJAC
*
      SHP7(4,1) = (- DYDETA*DY1 + DYDQSI*DY3) / D2JAC
      SHP7(4,2) = (- DYDETA*DY1 + DYDQSI*DY4) / D2JAC
      SHP7(4,3) = (- DYDETA*DY2 + DYDQSI*DY4) / D2JAC
      SHP7(4,4) = (- DYDETA*DY2 + DYDQSI*DY3) / D2JAC
      SHP7(5,1) = (  DYDETA*DX1 - DYDQSI*DX3) / D2JAC
      SHP7(5,2) = (  DYDETA*DX1 - DYDQSI*DX4) / D2JAC
      SHP7(5,3) = (  DYDETA*DX2 - DYDQSI*DX4) / D2JAC
      SHP7(5,4) = (  DYDETA*DX2 - DYDQSI*DX3) / D2JAC
*
      SHP8(4,1) = (  DXDETA*DY1 - DXDQSI*DY3) / D2JAC
      SHP8(4,2) = (  DXDETA*DY1 - DXDQSI*DY4) / D2JAC
      SHP8(4,3) = (  DXDETA*DY2 - DXDQSI*DY4) / D2JAC
      SHP8(4,4) = (  DXDETA*DY2 - DXDQSI*DY3) / D2JAC
      SHP8(5,1) = (- DXDETA*DX1 + DXDQSI*DX3) / D2JAC
      SHP8(5,2) = (- DXDETA*DX1 + DXDQSI*DX4) / D2JAC
      SHP8(5,3) = (- DXDETA*DX2 + DXDQSI*DX4) / D2JAC
      SHP8(5,4) = (- DXDETA*DX2 + DXDQSI*DX3) / D2JAC
*
      DO 300 I=1,NBNO
         K = 6*(I-1)
         IF (NOPLAN .EQ. 0) THEN
           BGR(1,K+1)=SHP(2,I)
           BGR(5,K+2)=SHP(3,I)
           BGR(1,K+5) = SHP(2,I)  * EXCEN
           BGR(5,K+4) = -SHP(3,I) * EXCEN
         ENDIF
         BGR(7,K+3)=SHP(2,I)
         BGR(3,K+4) = SHP7(4,I)
         BGR(3,K+5) = SHP7(5,I)
         BGR(8,K+3)=SHP(3,I)
         BGR(6,K+4) = SHP8(4,I)
         BGR(6,K+5) = SHP8(5,I)
*
         IF (.NOT. (IGAU .LE. 4)) THEN
            BGR(1,K+1)=SHP(2,I)
            BGR(5,K+2)=SHP(3,I)
            BGR(1,K+5) = SHP(2,I)  * EXCEN
            BGR(5,K+4) = -SHP(3,I) * EXCEN
            BGR(2,K+1)=SHP(3,I)
            BGR(4,K+2)=SHP(2,I)
            BGR(4,K+4) = -SHP(2,I) * EXCEN
            BGR(2,K+5) = SHP(3,I)  * EXCEN
         ENDIF
  300 CONTINUE
*
* On finit de completer la matrice BGR, pour les contraintes ou
* les deformations, par des valeurs calculees au centre:
*
      IF (IGAU .LE. 4) THEN
         DO 500 I=1,NBNO
            SHP(1,I)=SHPTOT(1,I, 5  )
            SHP(2,I)=SHPTOT(2,I, 5  )
            SHP(3,I)=SHPTOT(3,I, 5  )
  500    CONTINUE
*
* FONCTIONS DE FORME DANS LA GEOMETRIE REELLE:
*
         CALL JACOBI(XEL,SHP,2,NBNO,DJAC5)
         IF(DJAC5 .LE. ZER) THEN
            IARR=1
            RETURN
         END IF
*
         DO 520 I=1,NBNO
            K = 6*(I-1)
            IF (NOPLAN .EQ. 1) THEN
               BGR(1,K+1)=SHP(2,I)
               BGR(5,K+2)=SHP(3,I)
               BGR(1,K+5) = SHP(2,I)  * EXCEN
               BGR(5,K+4) = -SHP(3,I) * EXCEN
            ENDIF
            BGR(2,K+1)=SHP(3,I)
            BGR(4,K+2)=SHP(2,I)
            BGR(4,K+4) = -SHP(2,I) * EXCEN
            BGR(2,K+5) = SHP(3,I)  * EXCEN
  520    CONTINUE
      ENDIF
      RETURN
      END