Sources Esope | Cast3M

jacobix
C JACOBIX   SOURCE    AS218369  10/01/29    21:15:07     6623
C=======================================================================
C=                          J A C O B I X                              =
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=   TABA0(IDIM,NBNO) (E) Valeurs des ddl de saut aux noeuds           =
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 JACOBIX (XEL,TABA0,TABDH,SHP,IDIM,NBNO,DJac)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
      DIMENSION XEL(3,*),SHP(6,*)
      DIMENSION TABA0(IDIM,*),TABDH(*)
 
      PARAMETER (XZero=0.d0,XUn=1.d0)
 
      dJInv=XZero
 
C ===========================
C  1 - Cas de la DIMENSION 1
C ===========================
      IF (IDIM.EQ.1) THEN
        write(6,*) 'On ne traite pas la dimension 1 en XFEM'
        return
*        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))
          dXdQsi=dXdQsi+SHP(2,i)*(XEL(1,i)+TABDH(i)*TABA0(1,i))
          dXdEta=dXdEta+SHP(3,i)*(XEL(1,i)+TABDH(i)*TABA0(1,i))
          dYdQsi=dYdQsi+SHP(2,i)*(XEL(2,i)+TABDH(i)*TABA0(2,i))
          dYdEta=dYdEta+SHP(3,i)*(XEL(2,i)+TABDH(i)*TABA0(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
* 14/01/10: Il faut encore tester qu'on obtient de bons résultats !!!
        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)+TABDH(i)*TABA0(1,i))
          D21=D21+SHP(3,i)*(XEL(1,i)+TABDH(i)*TABA0(1,i))
          D31=D31+SHP(4,i)*(XEL(1,i)+TABDH(i)*TABA0(1,i))
          D12=D12+SHP(2,i)*(XEL(2,i)+TABDH(i)*TABA0(2,i))
          D22=D22+SHP(3,i)*(XEL(2,i)+TABDH(i)*TABA0(2,i))
          D32=D32+SHP(4,i)*(XEL(2,i)+TABDH(i)*TABA0(2,i))
          D13=D13+SHP(2,i)*(XEL(3,i)+TABDH(i)*TABA0(3,i))
          D23=D23+SHP(3,i)*(XEL(3,i)+TABDH(i)*TABA0(3,i))
          D33=D33+SHP(4,i)*(XEL(3,i)+TABDH(i)*TABA0(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
      ENDIF
 
      RETURN
      END