Esope source files

   1 : C CALJ22    SOURCE    CHAT      05/01/12    21:46:34     5004           
   2 :       SUBROUTINE CALJ22(FN,GR,PG,XYZ,HR,PGSQ,RPG,
   3 :      *ND,NP,NPG,IAXI,AIRE,AJ)
   4 : C************************************************************************
   5 : C
   6 : C          CE SP DIFFERRE DE CALI22 PAR LE RANGEMENT DE HR
   7 : C               HR(ND,NP)
   8 : C
   9 : C    REMARQUE IMPORTANTE : CE SP NE TRAITE QUE LES ELEMENTS VOLUMIQUES
  10 : C                           => HR(ND ...
  11 : C
  12 : C
  13 : C     CALCUL DE L'INVERSE DU JACOBIEN  AJ=1/J
  14 : C     CALCUL DE L'AIRE OU VOLUME       AIRE
  15 : C     CALCUL DE                        PGSQ(L)
  16 : C     CALCUL DE                        RPG(L)
  17 : C     CALCUL DE DES GRADIENTS          HR(ND,NP)
  18 : C     CALCUL INTERMEDIAIRE DE L'ELEMENT D'AIRE       SQ=DET(J)
  19 : C     DANS LES CAS 2D ET 3D
  20 : C
  21 : C     ND   DIMENSION ESPACE
  22 : C     NP   NOMBRE DE NOEUDS DE L'ELEMENT
  23 : C     NPG  NOMBRE DE POINTS D'INTEGRATION
  24 : C
  25 : C     XYZ  COORDONNEES
  26 : C     GR   GRADIENT
  27 : C     B & KAUX  TABLEAUX DE TRAVAIL
  28 : C************************************************************************
  29 :       IMPLICIT INTEGER(I-N)
  30 :       IMPLICIT REAL*8 (A-H,O-Z)
  31 : C
  32 :       REAL*8 FN(NP,NPG),GR(ND,NP,NPG),HR(ND,NP,NPG)
  33 :       REAL*8 PG(NPG),XYZ(ND,NP),PGSQ(NPG),RPG(NPG)
  34 :       REAL*8 SQ(64),B(3),AJ(ND,ND,NPG)
  35 :       DIMENSION KAUX(3)
  36 : -INC CCREEL
  37 : C
  38 : C***
  39 : 
  40 :       DEUPI=2.D0*XPI
  41 : 
  42 : C     WRITE(6,*)' SUB CALJ22 XYZ=',((XYZ(MM1,MM2),MM1=1,ND),MM2=1,NP)
  43 : 
  44 :       DO 10 L=1,NPG
  45 :       DO 10 M=1,ND
  46 :       DO 10 N=1,ND
  47 :       AJT=0.D0
  48 :       DO 11 I=1,NP
  49 :       AJT=AJT+GR(M,I,L)*XYZ(N,I)
  50 :    11 CONTINUE
  51 :       AJ(N,M,L)=AJT
  52 :    10 CONTINUE
  53 : C
  54 :       DO 20 L=1,NPG
  55 :       IF(ND.EQ.1)THEN
  56 :       VINT=AJ(1,1,L)
  57 : C     VINT=ABS(VINT)
  58 :       AJ(1,1,L)=1.D0/VINT
  59 :       ELSEIF(ND.EQ.2)THEN
  60 :       VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L)
  61 : C     VINT=ABS(VINT)
  62 :       RVINT=1.D0/VINT
  63 :       D11=AJ(2,2,L)
  64 :       D12=AJ(1,2,L)
  65 :       D21=AJ(2,1,L)
  66 :       D22=AJ(1,1,L)
  67 :       AJ(1,1,L)= RVINT*D11
  68 :       AJ(1,2,L)=-RVINT*D12
  69 :       AJ(2,1,L)=-RVINT*D21
  70 :       AJ(2,2,L)= RVINT*D22
  71 :       ELSEIF(ND.EQ.3)THEN
  72 :       D11=AJ(2,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(2,3,L)
  73 :       D12=AJ(1,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(1,3,L)
  74 :       D13=AJ(1,2,L)*AJ(2,3,L)-AJ(2,2,L)*AJ(1,3,L)
  75 :       D21=AJ(2,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(2,3,L)
  76 :       D22=AJ(1,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(1,3,L)
  77 :       D23=AJ(1,1,L)*AJ(2,3,L)-AJ(2,1,L)*AJ(1,3,L)
  78 :       D31=AJ(2,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(2,2,L)
  79 :       D32=AJ(1,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(1,2,L)
  80 :       D33=AJ(1,1,L)*AJ(2,2,L)-AJ(2,1,L)*AJ(1,2,L)
  81 :       VINT=AJ(1,1,L)*D11-AJ(1,2,L)*D21+AJ(1,3,L)*D31
  82 : C     VINT=ABS(VINT)
  83 :       RVINT=1.D0/VINT
  84 :       AJ(1,1,L)= RVINT*D11
  85 :       AJ(1,2,L)=-RVINT*D12
  86 :       AJ(1,3,L)= RVINT*D13
  87 :       AJ(2,1,L)=-RVINT*D21
  88 :       AJ(2,2,L)= RVINT*D22
  89 :       AJ(2,3,L)=-RVINT*D23
  90 :       AJ(3,1,L)= RVINT*D31
  91 :       AJ(3,2,L)=-RVINT*D32
  92 :       AJ(3,3,L)= RVINT*D33
  93 :       ELSE
  94 :       VINT=0.D0
  95 :       ENDIF
  96 : 
  97 :       SQ(L)=ABS(VINT)
  97 : C***  CALL DLAIN(ND,ND,AJ(1,1,L),KAUX,B,IER)
  98 : C
  99 :    20 CONTINUE
 100 : C
 101 :       DO 31 L=1,NPG
 102 :       DO 31 I=1,NP
 103 :       DO 31 N=1,ND
 104 :       U=0.D0
 105 :       DO 32 M=1,ND
 106 :    32 U=U+AJ(M,N,L)*GR(M,I,L)
 107 :    31 HR(N,I,L)=U
 108 :       U=0.D0
 109 :       DO 4 L=1,NPG
 110 :       U=U+PG(L)*SQ(L)
 111 :     4 PGSQ(L)=PG(L)*SQ(L)
 112 :       AIRE=U
 113 :       IF(IAXI.EQ.0)RETURN
 114 : C
 115 :       ID=3-IAXI
 115 :       IF(IAXI.EQ.3)CALL ARRET(0)
 116 :       DO 45 L=1,NPG
 117 :       RPGT  =0.D0
 118 :       DO 46 I=1,NP
 119 :       RPGT=RPGT+XYZ(ID,I)*FN(I,L)
 120 :  46   CONTINUE
 121 :       RPG(L)=RPGT
 122 :  45   CONTINUE
 123 : C
 124 :       U=0.D0
 125 :       DO 47 L=1,NPG
 126 :       PGSQ(L)=PGSQ(L)*DEUPI*RPG(L)
 127 :       U=U+PGSQ(L)
 128 :  47   CONTINUE
 129 :       AIRE=U
 130 : C
 131 :       RETURN
 132 :  1002 FORMAT(10(1X,1PE11.4))
 133 :  1001 FORMAT(20(1X,I5))
 134 :       END
 135 : 
 136 : 
 137 : 
 138 :  
 139 :  
 140 :