C JACOBI    SOURCE    CHAT      05/01/13    00:48:25     5004           

C=======================================================================
C=                          J A C O B I                                =
C=                          -----------                                =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=   Calcul du jacobien (DJac) et des derivees des fonctions de forme  =
C=   par rapport aux coordonnees REELLES x,y,z (SHPTOT(2 a 4,*)) en un =
C=   point donne d'un element (en general, un point de Gauss)          =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   XEL(3,NBNO)  (E)  Coordonnees GLOBALES des noeuds de l element    =
C=   SHP(6,NBNO) (E/S) Valeurs des fonctions de forme et de leurs      =
C=                     derivees au point considere (point de Gauss)    =
C=   IDIM         (E)  Dimension du probleme traite (1 a 3)            =
C=                     Cas particulier de l'element JOI3 : IDIM = 86   =
C=   NBNO         (E)  Nombre de NOEUDS de l element fini              =
C=   DJac         (S)  Jacobien calcule en ce point de Gauss           =
C=                     (egal a 0 si probleme en ce point)              =
C=                                                                     =
C=  Remarque :                                                         =
C=  ----------                                                         =
C=   Lors de l'entree dans le sous-programme, SHP(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 JACOBI (XEL,SHP,IDIM,NBNO,DJac)

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

      DIMENSION XEL(3,*),SHP(6,*)

      PARAMETER (XZero=0.,XUn=1.)

      dJInv=XZero

C ===========================
C  1 - Cas de la DIMENSION 1
C ===========================
      IF (IDIM.EQ.1) THEN
        DJac=XZero
        DO i=1,NBNO
          DJac=DJac+SHP(2,i)*XEL(1,i)
        ENDDO
        IF (DJac.NE.XZero) dJInv=XUn/DJac
        DO i=1,NBNO
          SHP(2,i)=SHP(2,i)*dJInv
        ENDDO
C ===========================
C  2 - Cas de la DIMENSION 2
C ===========================
      ELSE IF (IDIM.EQ.2) THEN
        dXdQsi=XZero
        dXdEta=XZero
        dYdQsi=XZero
        dYdEta=XZero
        DO 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)
        ENDDO
        DJac=dXdQsi*dYdEta-dXdEta*dYdQsi
        IF (DJac.NE.XZero) dJInv=XUn/DJac
        dQsidX= dYdEta*dJInv
        dQsidY=-dXdEta*dJInv
        dEtadX=-dYdQsi*dJInv
        dEtaDY= dXdQsi*dJInv
        DO i=1,NBNO
          V1=SHP(2,i)*dQsidX+SHP(3,i)*dEtadX
          SHP(3,i)=SHP(2,i)*dQsidY+SHP(3,i)*dEtadY
          SHP(2,i)=V1
        ENDDO
C ===========================
C  3 - Cas de la DIMENSION 3
C ===========================
      ELSE IF (IDIM.EQ.3) THEN
        D11=XZero
        D21=XZero
        D31=XZero
        D12=XZero
        D22=XZero
        D32=XZero
        D13=XZero
        D23=XZero
        D33=XZero
        DO i=1,NBNO
          D11=D11+SHP(2,i)*XEL(1,i)
          D21=D21+SHP(3,i)*XEL(1,i)
          D31=D31+SHP(4,i)*XEL(1,i)
          D12=D12+SHP(2,i)*XEL(2,i)
          D22=D22+SHP(3,i)*XEL(2,i)
          D32=D32+SHP(4,i)*XEL(2,i)
          D13=D13+SHP(2,i)*XEL(3,i)
          D23=D23+SHP(3,i)*XEL(3,i)
          D33=D33+SHP(4,i)*XEL(3,i)
        ENDDO
        DInv11=D22*D33-D23*D32
        DInv12=D32*D13-D12*D33
        DInv13=D12*D23-D22*D13
        DJac=D11*DInv11+D21*DInv12+D31*DInv13
        IF (DJac.NE.XZero) dJInv=XUn/DJac
        DInv11=DInv11*dJInv
        DInv12=DInv12*dJInv
        DInv13=DInv13*dJInv
        DInv21=(D23*D31-D21*D33)*dJInv
        DInv22=(D11*D33-D13*D31)*dJInv
        DInv23=(D21*D13-D11*D23)*dJInv
        DInv31=(D21*D32-D22*D31)*dJInv
        DInv32=(D12*D31-D11*D32)*dJInv
        DInv33=(D11*D22-D12*D21)*dJInv
        DO i=1,NBNO
          V1=DInv11*SHP(2,i)+DInv12*SHP(3,i)+DInv13*SHP(4,i)
          V2=DInv21*SHP(2,i)+DInv22*SHP(3,i)+DInv23*SHP(4,i)
          V3=DInv31*SHP(2,i)+DInv32*SHP(3,i)+DInv33*SHP(4,i)
          SHP(2,i)=V1
          SHP(3,i)=V2
          SHP(4,i)=V3
        ENDDO
C =======================================
C  4 - Cas particulier de l'element JOI3
C =======================================
      ELSE IF (IDIM.EQ.86) THEN
        dXdQsi=XZero
        dYdQsi=XZero
        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)
      ENDIF

      RETURN
      END