C BZCO      SOURCE    FANDEUR   13/01/16    21:15:00     7666
       SUBROUTINE BZCO(IGAU,MFR,IFOU,NIFOU,XEL,BPSS,SHPTOT,
     .                NST,NBNO,LRE,MELE,SHP,BGENE,DJAC,IRRT)
C=======================================================================
C     CALCUL DE LA MATRICE B
C     ET DU JACOBIEN EN IGAU POUR UN ELEMENT DE ZONE COHESIVE
C=======================================================================
C  INPUT
C     IGAU  = NUMERO DU POINT DE GAUSS
C     MFR   = NUMERO DE LA FORMULATION
C     XEL   = COORDONNEES DES NOEUDS DE L'ELEMENT
C     BPSS  = MATRICE DE PASSAGE
C     SHPTOT= FONCTIONS DE FORME ET DERIVEES DANS L'ESPACE DE REFERENCE
C             SHPTOT(1, ) = FONCTION DE FORME
C             SHPTOT(2, ) = DERIVEES PAR RAPPORT A QSI
C             SHPTOT(3, ) = DERIVEES PAR RAPPORT A ETA
C     NBNO  = NOMBRE DE NOEUDS
C     LRE   = NOMBRE DE COLONNES DE LA MATRICE B
C     NST   = NOMBRE DE COMPOSANTES DE CONTRAINTES
C  OUTPUT
C     SHP   = FONCTIONS DE FORME ET DERIVEES DANS L'ESPACE GEOMETRIQUE
C             SHP(1, ) = FONCTION DE FORME
C             SHP(2, ) = DERIVEES PAR RAPPORT A X LOCAL
C             SHP(3, ) = DERIVEES PAR RAPPORT A Y LOCAL
C     DJAC  = JACOBIEN AU POINT D'INTEGRATION IGAU
C     BGENE = MATRICE B AU POINT D'INTEGRATION IGAU
C     IRRT  = DIFFERENT DE ZERO SI ERREUR
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION XEL(3,*),BGENE(NST,*),SHP(6,*),SHPTOT(6,NBNO,*)
      DIMENSION BPSS(3,3)
C
      IRRT = 0
      NDIM = 2
      IF(IFOU .GT. 0) NDIM = 3
      IDIM = NDIM - 1
C
C     MATRICE JACOBIENNE
C
      DO 1 I=1,NBNO
       SHP(1,I) = SHPTOT(1,I,IGAU)
       SHP(2,I) = SHPTOT(2,I,IGAU)
    1 CONTINUE
C
C     TRAITEMENT PARTICULIER POUR LE CAS 2D
C     SINON, APPEL A DEVOLU
C
      IF(NDIM.EQ.2) THEN
        dXdQsi=0.D0
        dYdQsi=0.D0
        DO i=1,NBNO
          dXdQsi=dXdQsi+SHP(2,i)*XEL(1,i)
          dYdQsi=dYdQsi+SHP(2,i)*XEL(2,i)
        ENDDO
        DJAC=SQRT(dXdQsi*dXdQsi+dYdQsi*dYdQsi)
C
C  ON MULTIPLIE PAR LE RAYON EN AXI
C
       IF (IFOU .EQ.0) THEN
        RAYON=0.D0
        DO IRAY=1,NBNO
         RAYON=RAYON+SHPTOT(1,IRAY,IGAU)*XEL(1,IRAY)
        ENDDO
        DJAC=DJAC*RAYON
       ENDIF

      ELSE
        CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NIFOU,IDIM,1.D0,RR,DJAC)
      ENDIF

      IF (DJAC.LT.0.0D0) THEN
         IRRT = 1
      ELSE IF (DJAC.EQ.0.0D0) THEN
         IRRT = 2
      ENDIF
C
C     MATRICE B
C
      CALL ZERO(BGENE,NST,LRE)
C
      DO 2 I=1,NST
       DO 3 J=1,NBNO
        r_z = 2.D0 * SHP(1,J)
        L = NDIM*(J-1)
        DO 4 K=1,NDIM
         BGENE(I,L+K) = BPSS(I,K) * r_z
    4   CONTINUE
    3  CONTINUE
    2 CONTINUE
C
      RETURN
      END