C BELEC     SOURCE    FANDEUR   10/12/17    21:15:32     6427

C=======================================================================
C=                          B E L E C                                  =
C=                          ---------                                  =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=   Calcul de la matrice Bel reliant le champ electrique en un point  =
C=   d'un element fini aux ddls de potentiel aux noeuds de cet element =
C=   Le jacobien est egalement evalue en ce point pour verifier par la =
C=   suite si l'element fini n'est pas trop distordu.                  =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   XEL     (E)  Coordonnees des noeuds de l'element fini considere   =
C=   SHPREF  (E)  Fonctions de forme et leurs derivees au point etudie =
C=   NBNO    (E)  Nombre de noeuds de l'element fini                   =
C=   NDIM    (E)  Nombre de lignes de la matrice GRDEL                 =
C=   ISGNG   (E)  Valeur negative (<0) si on veut le champ electrique, =
C=                i.e. l'oppose du gradient                            =
C=                                                                     =
C=   SHP     (S)  Fonctions de forme et leurs derivees "actuelles"     =
C=   GRDEL   (S)  Matrice de gradients Bel calculee au point etudie    =
C=   DJAC    (S)  Jacobien au point de l'element etudie                =
C=                                                                     =
C=  Remarque :                                                         =
C=  ----------                                                         =
C=   Lors de l'entree dans le sous-programme, SHPREF(2 a 4,*) contient =
C=   les DERIVEES des fonctions de forme par rapport aux coordonnees   =
C=   de REFERENCE Qsi,Eta,Dzeta.                                       =
C=   En sortie du sous-programme, SHP(2 a 4,*) contient les DERIVEES   =
C=   des fonctions de FORME par rapport aux coordonnees REELLES x,y,z. =
C=======================================================================

      SUBROUTINE BELEC (XEL,SHPREF,NBNO,NDIM,ISGNG, SHP,GRDEL,DJAC)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC CCREEL

-INC PPARAM
-INC CCOPTIO

C= Quelques constantes (2.Pi et 4.Pi)
      PARAMETER (X2Pi= 6.283185307179586476925286766559D0)
      PARAMETER (X4Pi=12.566370614359172953850573533118D0)

      DIMENSION XEL(3,*),SHPREF(6,*), SHP(6,*),GRDEL(NDIM,*)

C ====================
C  1 - Initialisation
C ====================
C*OF     CALL ZERO(GRDEL,IDIM,NBNO)

C ===================================
C  2.1 - Cas des elements massifs 3D
C ===================================
      IF (IFOMOD.EQ.2) THEN
        DO j = 1, NBNO
          SHP(1,j) = SHPREF(1,j)
          SHP(2,j) = SHPREF(2,j)
          SHP(3,j) = SHPREF(3,j)
          SHP(4,j) = SHPREF(4,j)
        ENDDO
        CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
        IF (ISGNG .GE. 0) THEN
          DO j = 1, NBNO
            GRDEL(1,j) = SHP(2,j)
            GRDEL(2,j) = SHP(3,j)
            GRDEL(3,j) = SHP(4,j)
          ENDDO
        ELSE
          DO j = 1, NBNO
            GRDEL(1,j) = -SHP(2,j)
            GRDEL(2,j) = -SHP(3,j)
            GRDEL(3,j) = -SHP(4,j)
          ENDDO
        ENDIF
C =========================================================
C  2.2 - Cas des elements massifs 2D PLAN ou AXISymetrique
C =========================================================
      ELSE IF (IFOMOD.EQ.-1 .OR. IFOMOD.EQ.0) THEN
        DO j = 1,NBNO
          SHP(1,j) = SHPREF(1,j)
          SHP(2,j) = SHPREF(2,j)
          SHP(3,j) = SHPREF(3,j)
        ENDDO
        CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
        IF (ISGNG .GE. 0) THEN
          DO j = 1, NBNO
            GRDEL(1,j) = SHP(2,j)
            GRDEL(2,j) = SHP(3,j)
          ENDDO
        ELSE
          DO j = 1, NBNO
            GRDEL(1,j) = -SHP(2,j)
            GRDEL(2,j) = -SHP(3,j)
          ENDDO
        ENDIF
        IF (IFOMOD.EQ.0) THEN
          CALL DISTRR(XEL,SHP,NBNO,RR)
          DJAC = X2Pi*DJAC*RR
        ENDIF
C ===================================================
C  2.3 - Cas des elements massifs 2D en mode FOURIER
C ===================================================
      ELSE IF (IFOMOD.EQ.1) THEN
        DO j = 1,NBNO
          SHP(1,j) = SHPREF(1,j)
          SHP(2,j) = SHPREF(2,j)
          SHP(3,j) = SHPREF(3,j)
        ENDDO
        CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
        CALL DISTRR(XEL,SHP,NBNO,RR)
        IF (ISGNG .GE. 0) THEN
          DNN = NIFOUR/RR
          DO j = 1, NBNO
            GRDEL(1,j) = SHP(2,j)
            GRDEL(2,j) = DNN*SHP(1,j)
            GRDEL(3,j) = SHP(3,j)
          ENDDO
        ELSE
          DNN = -NIFOUR/RR
          DO j = 1, NBNO
            GRDEL(1,j) = -SHP(2,j)
            GRDEL(2,j) =  DNN*SHP(1,j)
            GRDEL(3,j) = -SHP(3,j)
          ENDDO
        ENDIF
        IF (NIFOUR.EQ.0) THEN
          DJAC = X2Pi*DJAC*RR
        ELSE
          DJAC = XPi*DJAC*RR
        ENDIF
C ===================================
C  2.4 - Cas des elements MASSIFS 1D
C ===================================
      ELSE IF (IFOMOD.EQ.3 .OR. IFOMOD.EQ.4 .OR. IFOMOD.EQ.5) THEN
        DO j = 1,NBNO
          SHP(1,j) = SHPREF(1,j)
          SHP(2,j) = SHPREF(2,j)
        ENDDO
        CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
C= Traitement des modes axisymetriques et spherique 1D
        IF (ISGNG .GE. 0) THEN
          DO j = 1,NBNO
            GRDEL(1,j) = SHP(2,j)
          ENDDO
        ELSE
          DO j = 1, NBNO
            GRDEL(1,j) = -SHP(2,j)
          ENDDO
        ENDIF
        IF (IFOMOD.EQ.4) THEN
          CALL DISTRR(XEL,SHP,NBNO,RR)
          DJAC = X2Pi*DJAC*RR
        ELSE IF (IFOMOD.EQ.5) THEN
          CALL DISTRR(XEL,SHP,NBNO,RR)
          DJAC = X4Pi*DJAC*RR*RR
        ENDIF

      ENDIF

      RETURN
      END