Sources Esope | Cast3M

caljtr
C CALJTR    SOURCE    CHAT      05/01/12    21:46:59     5004
      SUBROUTINE CALJTR(GR,XYZ,NES,ND,NP,NPG,AJ)
C************************************************************************
C
C     CALCUL DE L'INVERSE DU JACOBIEN TRANSPOSE QUI SERVIRA A CALCULER
C           t -1
C     û = (J )   u
C
C     CE SP TRAITE LES ELEMENTS VOLUMIQUES, SURFACIQUES et FILAIRES
C     DANS LES CAS 2D ET 3D
C
C     CALCUL TANGENTES ET NORMALES POUR LES ELEMENTS SURFACIQUES et
C            FILAIRES    : AJ=|tx ty|     |tx ty tz|
C                             |nx ny|  ou |ux uy uz|
C                                         |nx ny nz|
C
C
C     ND   DIMENSION ESPACE
C     NES  DIMENSION ESPACE ELEMENT DE REFERENCE
C     NP   NOMBRE DE NOEUDS DE L'ELEMENT
C     NPG  NOMBRE DE POINTS D'INTEGRATION
C
C     XYZ  COORDONNEES
C     GR   GRADIENT DANS L'ELEMENT DE REFERENCE
C
C************************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
      REAL*8 GR(NES,NP,NPG),XYZ(ND,NP),AJ(ND,ND,NPG)
C
C
C
C     write(6,*)' Debut CALJTR nes,nd,np,npg=',nes,nd,np,npg
      IF(NES.EQ.1.AND.ND.EQ.2)THEN
C
      DO 110 L=1,NPG
      AJX=0.D0
      AJY=0.D0
      DO 111 I=1,NP
      AJX=AJX+GR(1,I,L)*XYZ(1,I)
      AJY=AJY+GR(1,I,L)*XYZ(2,I)
  111 CONTINUE
      AJN=(AJX*AJX+AJY*AJY)**0.5D0
      AJ(1,1,L)=AJX/AJN
      AJ(2,1,L)=AJY/AJN
      AJ(1,2,L)=-AJY/AJN
      AJ(2,2,L)=AJX/AJN
 110  CONTINUE
C
      ELSEIF(NES.EQ.2.AND.ND.EQ.3)THEN
 
      AIRE=0.D0
      DO 210 L=1,NPG
      AJX=0.D0
      AJY=0.D0
      AJZ=0.D0
      BJX=0.D0
      BJY=0.D0
      BJZ=0.D0
      DO 211 I=1,NP
      AJX=AJX+GR(1,I,L)*XYZ(1,I)
      AJY=AJY+GR(1,I,L)*XYZ(2,I)
      AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
      BJX=BJX+GR(2,I,L)*XYZ(1,I)
      BJY=BJY+GR(2,I,L)*XYZ(2,I)
      BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
  211 CONTINUE
 
      XB=AJY*BJZ-AJZ*BJY
      YB=AJZ*BJX-AJX*BJZ
      ZB=AJX*BJY-AJY*BJX
 
      AJN=(XB*XB+YB*YB+ZB*ZB)**0.5D0
 
      AJ(1,1,L)=AJX/AJN
      AJ(2,1,L)=AJY/AJN
      AJ(3,1,L)=AJZ/AJN
      AJ(1,2,L)=BJX/AJN
      AJ(2,2,L)=BJY/AJN
      AJ(3,2,L)=BJZ/AJN
      AJ(1,3,L)=XB/AJN
      AJ(2,3,L)=YB/AJN
      AJ(3,3,L)=ZB/AJN
 210  CONTINUE
C
      ELSE
C
      DO 10 L=1,NPG
      DO 10 M=1,ND
      DO 10 N=1,ND
      AJT=0.D0
      DO 11 I=1,NP
      AJT=AJT+GR(M,I,L)*XYZ(N,I)
   11 CONTINUE
C     AJ(N,M,L)=AJT
C Jacobien Transposé
      AJ(M,N,L)=AJT
   10 CONTINUE
      ENDIF
C
      IF(ND.EQ.1)THEN
 
      DO 31 L=1,NPG
      VINT=AJ(1,1,L)
      AJ(1,1,L)=1.D0/VINT
 31   CONTINUE
 
      ELSEIF(ND.EQ.2)THEN
 
      DO 32 L=1,NPG
      VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L)
      RVINT=1.D0/VINT
      D11=AJ(2,2,L)
      D12=AJ(1,2,L)
      D21=AJ(2,1,L)
      D22=AJ(1,1,L)
      AJ(1,1,L)= RVINT*D11
      AJ(1,2,L)=-RVINT*D12
      AJ(2,1,L)=-RVINT*D21
      AJ(2,2,L)= RVINT*D22
 32   CONTINUE
 
      ELSEIF(ND.EQ.3)THEN
 
      DO 33 L=1,NPG
      D11=AJ(2,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(2,3,L)
      D12=AJ(1,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(1,3,L)
      D13=AJ(1,2,L)*AJ(2,3,L)-AJ(2,2,L)*AJ(1,3,L)
      D21=AJ(2,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(2,3,L)
      D22=AJ(1,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(1,3,L)
      D23=AJ(1,1,L)*AJ(2,3,L)-AJ(2,1,L)*AJ(1,3,L)
      D31=AJ(2,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(2,2,L)
      D32=AJ(1,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(1,2,L)
      D33=AJ(1,1,L)*AJ(2,2,L)-AJ(2,1,L)*AJ(1,2,L)
      VINT=AJ(1,1,L)*D11-AJ(1,2,L)*D21+AJ(1,3,L)*D31
      RVINT=1.D0/VINT
      AJ(1,1,L)= RVINT*D11
      AJ(1,2,L)=-RVINT*D12
      AJ(1,3,L)= RVINT*D13
      AJ(2,1,L)=-RVINT*D21
      AJ(2,2,L)= RVINT*D22
      AJ(2,3,L)=-RVINT*D23
      AJ(3,1,L)= RVINT*D31
      AJ(3,2,L)=-RVINT*D32
      AJ(3,3,L)= RVINT*D33
 33   CONTINUE
 
      ENDIF
C
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END