Sources Esope | Cast3M

jacobp
C JACOBP    SOURCE    CHAT      05/01/13    00:48:29     5004
      SUBROUTINE JACOBP(XEL,SHP,IDIM,NBBB,NBNO,DJAC)
*-----------------------------------------------------------------------
*
*     SOUS-PROGRAMME JACOBI ADAPTE AUX MILIEUX POREUX
*
*     CALCULE LES DERIVES DES FONCTIONS DE FORME
*     DANS LA GEOMETRIE DEFORMEE A PARTIR DES DERIVEES
*     DANS LA GEOMETRIE REDUITE AINSI QUE LE JACOBIEN
*   ENTREES
*     XEL(3,NBBB)=COORDONNEES LOCALES DE L ELEMENT
*     SHP(6,NBBB)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE
*     IDIM=DIMENSION
*     NBBB=NOMBRE DE NOEUDS
*     NBNO=NOMBRE DE FONCTIONS DE FORME
*   SORTIES
*     SHP(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMTRIE DE L ELEMENT
*     DJAC =JACOBIEN   ( MIS A 0 SI NON INVERSIBLE )
*
*-----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XEL(3,*),SHP(6,*)
      DIMENSION D(3,3),DINV(3,3),V(3)
      DXDQSI=0.D0
      DXDETA=0.D0
      DYDQSI=0.D0
      DYDETA=0.D0
C
C     CAS MONODIMENSIONEL
C
      IF(IDIM.NE.1) GOTO 20
      DO 400 I=1,NBBB
      DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
  400 CONTINUE
      DJAC=DXDQSI
      XXXX = DJAC
      IF(DJAC.NE.0.D0) XXXX=1.D0/DJAC
      DO 410 I=1,NBNO
      SHP(2,I)=SHP(2,I)*XXXX
  410 CONTINUE
      GOTO 666
C
C     CAS 2 DIMENSIONS
C
  20  CONTINUE
      IF (IDIM.NE.2) GOTO 30
      DO 100 I=1,NBBB
      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)
 100  CONTINUE
      DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
      XXXX = DJAC
      IF(DJAC.NE.0.D0) XXXX=1.D0/DJAC
      DQSIDX= DYDETA*XXXX
      DQSIDY=-DXDETA*XXXX
      DETADX=-DYDQSI*XXXX
      DETADY= DXDQSI*XXXX
      DO 200 I=1,NBNO
      TT=SHP(2,I)*DQSIDX+SHP(3,I)*DETADX
      SHP(3,I)=SHP(2,I)*DQSIDY+SHP(3,I)*DETADY
      SHP(2,I)=TT
  200 CONTINUE
      GOTO 666
C
C     CAS TRIDIMENSIONEL
C
  30  IF (IDIM.NE.3) GOTO 666
      DO 310 I=1,3
      DO 310 J=1,3
      D(I,J)=0.D0
 310  CONTINUE
      DO 300 I=1,NBBB
C
      D(1,1)=D(1,1)+SHP(2,I)*XEL(1,I)
      D(2,1)=D(2,1)+SHP(3,I)*XEL(1,I)
      D(3,1)=D(3,1)+SHP(4,I)*XEL(1,I)
C
      D(1,2)=D(1,2)+SHP(2,I)*XEL(2,I)
      D(2,2)=D(2,2)+SHP(3,I)*XEL(2,I)
      D(3,2)=D(3,2)+SHP(4,I)*XEL(2,I)
C
      D(1,3)=D(1,3)+SHP(2,I)*XEL(3,I)
      D(2,3)=D(2,3)+SHP(3,I)*XEL(3,I)
      D(3,3)=D(3,3)+SHP(4,I)*XEL(3,I)
 300  CONTINUE
C     INVERSION DE LA MATRICE D(3,3)
      DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2)
      DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1)
      DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1)
      DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3)
      DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1)
      DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3)
      DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1)
      DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3)
      DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1)
      DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3)
      XXXX = DJAC
      IF(DJAC.NE.0.D0) XXXX=1.D0/DJAC
      DO 350 IA=1,3
      DO 350 IB=1,3
      DINV(IA,IB)=DINV(IA,IB)*XXXX
  350 CONTINUE
      DO 360 IA=1,NBNO
      DO 370 IB=1,3
      V(IB)=0.D0
      DO 370 IC=1,3
      V(IB)=V(IB)+DINV(IB,IC)*SHP(IC+1,IA)
  370 CONTINUE
      SHP(2,IA)=V(1)
      SHP(3,IA)=V(2)
      SHP(4,IA)=V(3)
 360  CONTINUE
  666 CONTINUE
      RETURN
      END