Sources Esope | Cast3M

caljcc
C CALJCC    SOURCE    CHAT      05/01/12    21:46:46     5004
      SUBROUTINE CALJCC(FN,GR,PG,XYZ,HR,PGSQ,RPG,IES,ND,NP,NPG,IAXI,
     &AIRE,AJ,HHR)
C************************************************************************
C
C          CE SP DIFFERRE DE CALICC PAR LE RANGEMENT DE HR
C
C                        HR(IES,NP,1)
C                        HHR(ND,NP,1)
C
C     DANS LE CAS DES ELEMENTS COQUES LE CHANGEMENT DE REPERE SE FAIT EN
C     TEMPS 1/ DANS LE PLAN DE L'ELEMENT : HHR
C           2/ ROTATION 3D               : HR
C
C     CALCUL DE L'INVERSE DU JACOBIEN  AJ=1/J
C     CALCUL DE L'AIRE OU VOLUME       AIRE
C     CALCUL DE                        PGSQ(L)
C     CALCUL DE                        RPG(L)
C     CALCUL DE DES GRADIENTS          HR(ND,NP)
C     CALCUL INTERMEDIAIRE DE L'ELEMENT D'AIRE       SQ=DET(J)
C     DANS LES CAS 2D ET 3D
C
C     IES  DIMENSION ESPACE
C     ND   DIMENSION ESPACE DE L'ELEMENT
C     NP   NOMBRE DE NOEUDS DE L'ELEMENT
C     NPG  NOMBRE DE POINTS D'INTEGRATION
C
C     XYZ  COORDONNEES
C     GR   GRADIENT
C     B & KAUX  TABLEAUX DE TRAVAIL
C************************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
-INC CCREEL
C
      REAL*8 FN(NP,NPG),GR(ND,NP,1),HR(IES,NP,1),HHR(ND,NP,1)
      REAL*8 PG(NPG),XYZ(IES,NP),PGSQ(NPG),RPG(NPG),XY(2,9)
      REAL*8 SQ(9),B(3),AJ(IES,IES,NPG),AJJ(3,3,9)
C
      DIMENSION KAUX(3)
C
C***
C     WRITE(6,*) ' *** SUB CALJCC *** IES=',IES
 
      IF(IES.NE.2)GO TO 30
 
C                ----------   2D   ----------
C CAS DES ELEMENTS POUTRE 2D
 
      TX=XYZ(1,NP)-XYZ(1,1)
      TY=XYZ(2,NP)-XYZ(2,1)
 
      AIRE=TX*TX+TY*TY
      AIRE=SQRT(AIRE)
 
      TX=TX/AIRE
      TY=TY/AIRE
 
      PX=-TY
      PY=TX
 
C     WRITE(6,*)' TX,TY=',TX,TY
C     WRITE(6,*)' PX,PY=',PX,PY
C     WRITE(6,*)' AIRE =',AIRE,'  NPG=',NPG
C     WRITE(6,*)' PG   =',PG
C     WRITE(6,*)' GR   =',GR
C     CALL ARRET(0)
 
 
 
      DO 20 L=1,NPG
 
      AJ(1,1,L)=TX
      AJ(1,2,L)=TY
      AJ(2,1,L)=PX
      AJ(2,2,L)=PY
      PGSQ(L)=PG(L)*AIRE
 
      DO 20 I=1,NP
      HR(1,I,L)=AJ(1,1,L)*GR(1,I,L)/PGSQ(L)
      HR(2,I,L)=AJ(1,2,L)*GR(1,I,L)/PGSQ(L)
 20   CONTINUE
 
      IF(IAXI.EQ.0)RETURN
C
      ID=3-IAXI
      IF(IAXI.EQ.3)CALL ARRET(0)
      DO 25 L=1,NPG
      RPG(L)=XZERO
      DO 26 I=1,NP
      RPG(L)=RPG(L)+XYZ(ID,I)*FN(I,L)
 26   CONTINUE
 25   CONTINUE
C
      AIRE=XZERO
      DO 27 L=1,NPG
      PGSQ(L)=PGSQ(L)*2.0D0*XPI*RPG(L)
      AIRE=AIRE+PGSQ(L)
 27   CONTINUE
      RETURN
 
 30   CONTINUE
C                ----------   3D   ----------
      IF(ND.NE.1)GO TO 40
C CAS DES ELEMENTS POUTRE 3D
 
      TX=XYZ(1,NP)-XYZ(1,1)
      TY=XYZ(2,NP)-XYZ(2,1)
      TZ=XYZ(3,NP)-XYZ(3,1)
 
      AIRE=TX*TX+TY*TY+TZ*TZ
      AIRE=SQRT(AIRE)
 
      TX=TX/AIRE
      TY=TY/AIRE
      TZ=TZ/AIRE
 
      QX=3*TY-2*TZ
      QY=TZ-3*TX
      QZ=2*TX-TY
 
      QQ=QX*QX+QY*QY+QZ*QZ
      QQ=SQRT(QQ)
 
      QX=QX/QQ
      QY=QY/QQ
      QZ=QZ/QQ
 
      PX=QY*TZ-QZ*TY
      PY=QZ*TX-QX*TZ
      PZ=QX*TY-QY*TX
 
 
      DO 31 L=1,NPG
 
      AJ(1,1,L)=TX
      AJ(1,2,L)=TY
      AJ(1,3,L)=TZ
      AJ(2,1,L)=PX
      AJ(2,2,L)=PY
      AJ(2,3,L)=PZ
      AJ(3,1,L)=QX
      AJ(3,2,L)=QY
      AJ(3,3,L)=QZ
      PGSQ(L)=PG(L)*AIRE
 
      DO 31 I=1,NP
C     HR(1,I,L)=AJ(1,1,L)*GR(1,I,L)
      HR(1,I,L)=GR(1,I,L)/AIRE
C     HR(2,I,L)=AJ(1,2,L)*GR(1,I,L)
      HR(2,I,L)=XZERO
C     HR(3,I,L)=AJ(1,3,L)*GR(1,I,L)
      HR(3,I,L)=XZERO
 31   CONTINUE
      RETURN
 
 40   CONTINUE
C     WRITE(6,*)' *** CAS DES ELEMENTS COQUES 3D ***'
 
 
      CALL CALJQB(XYZ,AJJ,3,NP)
      CALL CALJXY(XYZ,AJJ,IES,XY,ND,NP)
 
 
 
C     WRITE(6,*)' SUB CALJCC : APPEL A CALJ22 '
C     WRITE(6,*)' ND,NP,NPG,IAXI=',ND,NP,NPG,IAXI
 
C     DO 441 L=1,NPG
C     WRITE(6,*)' L=',L
C     WRITE(6,1002) (GR(1,I,L),I=1,NP)
C     WRITE(6,1002) (GR(2,I,L),I=1,NP)
C441  CONTINUE
 
      CALL CALJ22(FN,GR,PG,XY,HHR,PGSQ,RPG,ND,NP,NPG,IAXI,AIRE,AJ)
 
      DO 41 L=1,NPG
      PGSQ(L)=PG(L)*AIRE
      DO 41 I=1,NP
      DO 41 N=1,IES
      HR(N,I,L)=XZERO
      DO 41 M=1,ND
      HR(N,I,L)=HR(N,I,L)+AJJ(M,N,1)*HHR(M,I,L)
 41   CONTINUE
 
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END