Sources Esope | Cast3M

jacobi
C JACOBI    SOURCE    PV090527  23/07/13    21:15:04     11708          
 
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 (XUn=1.d0)
-INC CCREEL
 
      dJInv=xgrand*xzprec
 
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
          shp2=shp(2,i)
          shp3=shp(3,i)
          shp4=shp(4,i)
          xel1=xel(1,i)
          xel2=xel(2,i)
          xel3=xel(3,i)
          D11=D11+SHP2*XEL1
          D21=D21+SHP3*XEL1
          D31=D31+SHP4*XEL1
          D12=D12+SHP2*XEL2
          D22=D22+SHP3*XEL2
          D32=D32+SHP4*XEL2
          D13=D13+SHP2*XEL3
          D23=D23+SHP3*XEL3
          D33=D33+SHP4*XEL3
        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
          shp2=shp(2,i)
          shp3=shp(3,i)
          shp4=shp(4,i)
          V1=DInv11*SHP2+DInv12*SHP3+DInv13*SHP4
          V2=DInv21*SHP2+DInv22*SHP3+DInv23*SHP4
          V3=DInv31*SHP2+DInv32*SHP3+DInv33*SHP4
          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