1 : C KSPRUS SOURCE MAGN 10/05/18 21:16:08 6675 2 : SUBROUTINE KSPRUS 3 : &(MELEME,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0, 4 : & ANUKK,IDCEN,IKOMP,NPTU,IPADU,IZTU1,TETAN,NP,NEL,CMD) 5 : 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8 (A-H,O-Z) 8 : C***************************************************************************** 9 : C 10 : C Ce SP calcule les matrices elementaires de divergence alias C 11 : C 12 : C IKAS=1 KMCT calcul de Ct (Div U) 13 : C IKAS=2 KMAC calcul de C uniquement (Grad p) 14 : C IKAS=3 KCCT calcul de C assemblage pour C et Ct 15 : C 16 : C***************************************************************************** 17 : CHARACTER*8 NOM0 18 : 19 : -INC CCOPTIO 20 : -INC CCGEOME 21 : -INC SMCOORD 22 : -INC SIZFFB 23 : POINTEUR IZF1.IZFFM 24 : -INC SMMATRIK 25 : -INC SMLENTI 26 : POINTEUR IPADU.MLENTI 27 : -INC SMELEME 28 : -INC SMCHPOI 29 : POINTEUR IZTGG1.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL 30 : -INC CCREEL 31 : PARAMETER (LRV=64) 32 : SEGMENT PETROV 33 : REAL*8 WT(LRV,NP,NPG,IDIM),WS(LRV,NP,NPG,IDIM),HK(LRV,IDIM,NP,NPG) 34 : REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) 35 : REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV) 36 : REAL*8 AJK(LRV,IDIM,IDIM,NPG) 37 : ENDSEGMENT 38 : DIMENSION KIPM(3),XYZ1(24) 39 : C 40 : C OPERATEUR PRESSION 41 : C 42 : DEUPI=1.D0 43 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 44 : IF(IKAS.EQ.3)DEUPI=-DEUPI 45 : 46 : NBEL = NEL 47 : C write(6,*)' DEUPI=',deupi 48 : IF(IDIM.EQ.2)IPM3=IPM1 49 : KIPM(1)=IPM1 50 : KIPM(2)=IPM2 51 : KIPM(3)=IPM3 52 : SEGACT MELEME 53 : SEGACT IPM1*MOD,IPM2*MOD,IPM3*MOD 54 : SEGACT IPADU,IZTU1,TETAN,IZTGG1 55 : 56 : IF(KPRE.NE.2)THEN 57 : C 58 : C ON TRAITE LE CAS DIV U -> T sommet 59 : C 60 : IF(INEFMD.EQ.3)THEN 61 : IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'PRP0' 62 : IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'PRP1' 63 : IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'PFP1' 64 : ELSEIF(INEFMD.EQ.2)THEN 65 : IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'MCP0' 66 : IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'MCP1' 67 : IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'MCF1' 68 : ELSEIF(INEFMD.EQ.4)THEN 69 : NOM0=NOMS(ITYPEL)//' ' 70 : ENDIF 71 : c write(*,*)' Kpruss : NOM0=',nom0,' ikas=',ikas,' nel=',nel 72 : c write(*,*)' IK1',IK1,' COEF',IZTGG1.VPOCHA(1,1),'ANUK',ANUKK 72 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 73 : 74 : SEGACT IZFFM*MOD 75 : IZHR=KZHR(1) 76 : SEGACT IZHR*MOD 77 : NES=GR(/1) 78 : NPG=GR(/3) 79 : IZF1=KTP(1) 80 : SEGACT IZF1*MOD 81 : * ON modifie pour avoir les bonnes fonctions de forme. GBM 18/10/99 82 : MP1=FN(/1) 83 : * 84 : * On fait des packets 85 : * 86 : NNN=MOD(NBEL,LRV) 87 : IF(NNN.EQ.0) THEN 88 : NPACK=NBEL/LRV 89 : ELSE 90 : NPACK=1+(NBEL-NNN)/LRV 91 : ENDIF 92 : KPACKD=1 93 : KPACKF=NPACK 94 : * 95 : * boucle sur les packets, on active la structure petrov 96 : * 97 : SEGINI PETROV 98 : DO 7001 KPACK=KPACKD,KPACKF 99 : C 100 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 101 : C 102 : C 1. Calcul des limites du paquet courant. 103 : KDEB=1+(KPACK-1)*LRV 104 : KFIN=MIN(NBEL,KDEB+LRV-1) 105 : C 106 : DO 7002 K=KDEB,KFIN 107 : KP=K-KDEB+1 108 : NK = K + K0 109 : NKR=(1-IK1)*(NK-1)+1 110 : COEFK(KP)=IZTGG1.VPOCHA(NKR,1) + 1.D-30 111 : * WRITE(*,*) IZTGG1.VPOCHA(NKR,1), COEFK(KP) 112 : ANUK(KP)=ANUKK/COEFK(KP) 113 : 7002 CONTINUE 114 : C 115 : C Une seule composante de projection car div se projete 116 : C sur un scalaire. 117 : C 118 : c write (*,*) 'IPADU', IPADU.LECT(/1) 118 : CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 119 : & NES,NP,NPG,IAXI,XCOOR, 120 : & MELEME.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, 121 : & TETAN.VPOCHA,IPADU.LECT,IZTU1.VPOCHA,IPADU.LECT,NPTU,ANUK, 122 : & WT(1,1,1,1),WS(1,1,1,1),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, 123 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 124 : 125 : DO 30 KE=KDEB,KFIN 126 : NK=KE-KDEB+1 127 : 128 : * DO 12 I=1,NP 129 : * J=NUM(I,KE) 130 : * DO 12 N=1,IDIM 131 : * XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N) 132 : * 12 CONTINUE 133 : 133 : * CALL CALJBR 134 : * &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 135 : C write(6,*)' Retour caljbr ',mp1,np,npg 136 : 137 : 138 : 139 : DO 324 K=1,IDIM 140 : IPM4=KIPM(K) 141 : DO 323 M=1,MP1 142 : DO 323 I=1,NP 143 : U=0.D0 144 : DO 333 L=1,NPG 145 : * ON modifie la formule GBM 18/10/99 146 : U=U+WT(NK,M,L,1)*HK(NK,K,I,L)*PGSK(NK,L)*DEUPI*RPGK(NK,L) 147 : * U=U+FN(M,L)*HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L) 148 : 333 CONTINUE 149 : 150 : IF(IAXI.NE.0.AND.K.EQ.1)THEN 151 : DO 334 L=1,NPG 152 : * modifié par GBM 18/10/99 153 : U=U+WT(NK,M,L,1)*FN(I,L)*PGSK(NK,L)*DEUPI 154 : * U=U+FN(M,L)*FN(I,L)*PGSQ(L)*DEUPI 155 : 334 CONTINUE 156 : ENDIF 157 : C 158 : C K1=1+(1-IK1)*(NK-1) 159 : C attention valable que pour coef au centre 160 : U=U*COEFK(NK) 161 : * WRITE(*,*) NK,M,I,COEFK(NK) 162 : C 163 : if(ikas.ne.2)then 164 : IPM4.AM(KE,I,M)=IPM4.AM(KE,I,M)+U 165 : else 166 : IPM4.AM(KE,M,I)=IPM4.AM(KE,M,I)+U 167 : endif 168 : 169 : * write(*,*) KE,PGSQ 170 : 323 CONTINUE 171 : 324 CONTINUE 172 : 30 CONTINUE 173 : 7001 CONTINUE 174 : 175 : SEGSUP IZHR,IZFFM,PETROV 176 : 177 : 178 : C CAS MACRO CENTRE 179 : 180 : ELSEIF(KPRE.EQ.2)THEN 181 : * initialisation de npg au pif par TC 182 : NPG=1 183 : SEGINI PETROV 184 : NOM0=NOMS(ITYPEL)//' ' 184 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 185 : 186 : SEGACT IZFFM*MOD 187 : IZHR=KZHR(1) 188 : SEGACT IZHR*MOD 189 : NES=GR(/1) 190 : NPG=GR(/3) 191 : IZF1=KTP(1) 192 : SEGACT IZF1*MOD 193 : MPG=IZF1.FN(/2) 194 : NP=GR(/2) 195 : 196 : DO 40 KE=1,NEL 197 : 198 : NK=K0+KE 199 : 200 : IX=0 201 : DO 42 I=1,NP 202 : J=MELEME.NUM(I,KE) 203 : DO 42 N=1,IDIM 204 : IX=IX+1 205 : XYZ1(IX)=XCOOR((J-1)*(IDIM+1) +N) 206 : 42 CONTINUE 207 : 207 : CALL CALJBR(FN,GR,PG,XYZ1,HR,PGSQ,RPG,NES, 208 : & IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 209 : 210 : DO 424 K=1,IDIM 211 : IPM4=KIPM(K) 212 : 213 : DO 423 I=1,NP 214 : U=0.D0 215 : DO 433 L=1,NPG 216 : U=U+HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L) 217 : 433 CONTINUE 218 : 219 : IF(IAXI.NE.0.AND.K.EQ.1)THEN 220 : DO 434 L=1,NPG 221 : U=U+FN(I,L)*PGSQ(L)*DEUPI 222 : 434 CONTINUE 223 : ENDIF 224 : 225 : K1=1+(1-IK1)*(NK-1) 226 : U=U*IZTGG1.VPOCHA(K1,1) 227 : 228 : if(ikas.ne.2)then 229 : IPM4.AM(KE,I,1)=IPM4.AM(KE,I,1)+U 230 : else 231 : IPM4.AM(KE,1,I)=IPM4.AM(KE,1,I)+U 232 : endif 233 : 234 : 423 CONTINUE 235 : 424 CONTINUE 236 : 237 : 40 CONTINUE 238 : 239 : SEGSUP IZHR,IZFFM,PETROV 240 : 241 : ENDIF 242 : 243 : RETURN 244 : 1002 FORMAT(10(1X,1PE11.4)) 245 : 1040 FORMAT(1X,'CALCUL MATRICE AM ',I4/10(1X,1PE11.4)) 246 : END 247 : 248 : 249 : 250 : 251 : 252 : 253 : 254 : 255 : 256 : 257 : 258 : 259 : 260 : 261 : 262 : 263 : 264 : 265 :
© Cast3M 2003 - All rights reserved.
Disclaimer