Sources Esope | Cast3M

calj22
C CALJ22    SOURCE    CHAT      05/01/12    21:46:34     5004
      SUBROUTINE CALJ22(FN,GR,PG,XYZ,HR,PGSQ,RPG,
     *ND,NP,NPG,IAXI,AIRE,AJ)
C************************************************************************
C
C          CE SP DIFFERRE DE CALI22 PAR LE RANGEMENT DE HR
C               HR(ND,NP)
C
C    REMARQUE IMPORTANTE : CE SP NE TRAITE QUE LES ELEMENTS VOLUMIQUES
C                           => HR(ND ...
C
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     ND   DIMENSION ESPACE
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
      REAL*8 FN(NP,NPG),GR(ND,NP,NPG),HR(ND,NP,NPG)
      REAL*8 PG(NPG),XYZ(ND,NP),PGSQ(NPG),RPG(NPG)
      REAL*8 SQ(64),B(3),AJ(ND,ND,NPG)
      DIMENSION KAUX(3)
-INC CCREEL
C
C***
 
      DEUPI=2.D0*XPI
 
C     WRITE(6,*)' SUB CALJ22 XYZ=',((XYZ(MM1,MM2),MM1=1,ND),MM2=1,NP)
 
      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
      AJ(N,M,L)=AJT
   10 CONTINUE
C
      DO 20 L=1,NPG
      IF(ND.EQ.1)THEN
      VINT=AJ(1,1,L)
C     VINT=ABS(VINT)
      AJ(1,1,L)=1.D0/VINT
      ELSEIF(ND.EQ.2)THEN
      VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L)
C     VINT=ABS(VINT)
      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
      ELSEIF(ND.EQ.3)THEN
      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
C     VINT=ABS(VINT)
      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
      ELSE
      VINT=0.D0
      ENDIF
 
      SQ(L)=ABS(VINT)
C***  CALL DLAIN(ND,ND,AJ(1,1,L),KAUX,B,IER)
C
   20 CONTINUE
C
      DO 31 L=1,NPG
      DO 31 I=1,NP
      DO 31 N=1,ND
      U=0.D0
      DO 32 M=1,ND
   32 U=U+AJ(M,N,L)*GR(M,I,L)
   31 HR(N,I,L)=U
      U=0.D0
      DO 4 L=1,NPG
      U=U+PG(L)*SQ(L)
    4 PGSQ(L)=PG(L)*SQ(L)
      AIRE=U
      IF(IAXI.EQ.0)RETURN
C
      ID=3-IAXI
      IF(IAXI.EQ.3)CALL ARRET(0)
      DO 45 L=1,NPG
      RPGT  =0.D0
      DO 46 I=1,NP
      RPGT=RPGT+XYZ(ID,I)*FN(I,L)
 46   CONTINUE
      RPG(L)=RPGT
 45   CONTINUE
C
      U=0.D0
      DO 47 L=1,NPG
      PGSQ(L)=PGSQ(L)*DEUPI*RPG(L)
      U=U+PGSQ(L)
 47   CONTINUE
      AIRE=U
C
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END