1 : C KRESS SOURCE PV 08/09/11 21:16:11 6150 2 : SUBROUTINE KRESS(MTABP,MCHB,MCHR,IMPR,BETA,KDPDQ,KPIMP,PIMP,NIMP) 3 : C *********************************************************************** 4 : C 5 : C RESOLUTION DE LA PRESSION DANS LE CAS D'UNE METHODE ITERATIVE 6 : C AVEC MATRICE STOCKEE. ON ENCHAINE, SI BESOIN EST : PROGCS, 7 : C CGPR ET RESCGS. 8 : C 9 : C DANS LE CAS D'UNE METHODE ITERATIVE AVEC STOCKAGE DE LA MATRICE 10 : C DE PRESSION .... 11 : C 12 : C CE SOUS PROGRAMME EFFECTUE LES OPERATIONS SUIVANTES : 13 : C 14 : C 1. STOCKAGE DANS LA TABLE METHODE DES POINTEURS ET DONNEES NECESSAIRES 15 : C 16 : C 2. CALCUL ET STOCKAGE DE LA MATRICE C D C CONNAISSANT C 17 : C EN FONCTION DU MODE DE STOCKAGE DEMANDE 18 : C 19 : C DEUX MODES DE STOCKAGE SONT DISPONIBLES 20 : C - MORSE CLASSIQUE (KSTOCK.EQ.0) 21 : C - MODE COMPRESSE (KSTOCK.EQ.1) 22 : C 23 : C 3. CALCUL D'UN PRECONDITIONNEMENT SI NECESSAIRE, I.E. LORSQUE 24 : C KTYPI=3 OU 4 25 : C 26 : C 27 : C 28 : C 29 : C POINTEURS : 30 : C 31 : C MTABP CONTIENT ENTRE AUTRES LES TABLEAUX D'INDIRECTION 32 : C (CF PROFIM) 33 : C MELEME OBJET MAILLAGE SUR LEQUEL EST FAIT LE CALCUL 34 : C MATRAK MATRICES ELEMENTAIRES DE LA DIVERGENCE (ALIAS"C") 35 : C 36 : C IZIPAD CORRESPONDANCE NUMER. GLOBALE --> NUMER. LOCALE 37 : C (DOMAINE SUR LEQUEL PORTE AP ET NON LA PRESSION) 38 : C IZTAB POROSITE (VECT DOMA) OU (VECT NOEU) PAR DEFAUT 1.,1.,1. 39 : C 40 : C EN SORTIE : 41 : C 42 : C LA MATRICE EST RANGEE DANS UN SEUL SEGMENT DE MEMOIRE 43 : C 44 : C MORSE : PAR LIGNES, LE PREMIER DE LA LIGNE EST L'ELEMENT LUI-MEME 45 : C COMPRESSE : COMME UN TABLEAU A(NEL,NNZ) AVEC NNZ = NOMBRE MAXI 46 : C DE TERMES NON NULS PAR LIGNE. LA DIAGONALE EST STOCKEE 47 : C EN PREMIER. 48 : C 49 : C LE CALCUL DE LA MATRICE EST FAIT SOUS-OBJET PAR SOUS-OBJET CE QUI 50 : C LIMITE A DEUX LE NOMBRE DE MELEMEs ET DE MATRICES DE LA DIVERGENCE 51 : C ACTIVES AU MEME INSTANT (EN PLUS DE LA MATRICE ET DE SON TABLEAU 52 : C D'INDEXAGE) 53 : C 54 : C *********************************************************************** 55 : IMPLICIT INTEGER(I-N) 56 : IMPLICIT REAL*8 (A-H,O-Z) 57 : C 58 : C*** 59 : 60 : -INC CCOPTIO 61 : -INC SMELEME 62 : POINTEUR IPTJ.MELEME,IPTK.MELEME,MELEM1.MELEME 63 : POINTEUR IGEOM.MELEME,IGEOMI.MELEME 64 : POINTEUR MELEMK.MELEME,MELEMC.MELEME,IGEOMC.MELEME 65 : POINTEUR MELSTB.MELEME 66 : -INC SMLENTI 67 : POINTEUR IZIPAD.MLENTI,MLENTP.MLENTI,IZIPAK.MLENTI 68 : POINTEUR IA.MLENTI,JA.MLENTI,KA.MLENTI 69 : -INC SMCOORD 70 : -INC SMCHPOI 71 : POINTEUR IZDV.MCHPOI 72 : POINTEUR IZDP.MCHPOI,IZDDP.MPOVAL 73 : POINTEUR MCHB.MCHPOI,MCHR.MCHPOI 74 : POINTEUR IZR.MPOVAL,IZPOC.MPOVAL 75 : POINTEUR IZP.MPOVAL,IZF.MPOVAL 76 : POINTEUR IPRESU.MCHPOI 77 : POINTEUR IPRESV.MPOVAL 78 : C-INC SMMATRAK 79 : -INC SMMATRK1 80 : C************************************************************************* 81 : C 82 : C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees 83 : C 84 : 85 : * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange 86 : * (points CENTRE ) pour chaque operateur de contrainte 87 : * KGEOC SPG pour la totalite des points CENTRE. 88 : * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse) 89 : * KLEMC Connectivites de l'ensemble des contraintes 90 : * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones 91 : 92 : SEGMENT MATRAK 93 : INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP) 94 : INTEGER LIZAFM(NBSOUS) 95 : INTEGER IKAM0 (NBSOUS) 96 : INTEGER IMEM (NBELC) 97 : INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC 98 : ENDSEGMENT 99 : 100 : SEGMENT IZAFM 101 : REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX) 102 : ENDSEGMENT 103 : 104 : POINTEUR IPMJ.IZAFM,IPMK.IZAFM 105 : 106 : C******************************************************************* 107 : -INC SMLREEL 108 : POINTEUR IZB.MLREEL 109 : POINTEUR IZDPR.MLREEL,IZPIM.MLREEL 110 : POINTEUR IZICCG.MLREEL,IZA.MLREEL 111 : POINTEUR IZD.MLREEL 112 : C 113 : CHARACTER*8 TYPE,TYPC 114 : C 115 : C Pour verifier ou non les numeros des elements sur 116 : C une condition de sortie fluide. 117 : C 118 : LOGICAL SOMTER 119 : DATA SOMTER/.FALSE./ 120 : 121 : C*** 122 : IZL=0 123 : KDPDQ=0 124 : C write(6,*)' Entree dans kress KDPDQ=',kdpdq 125 : 126 : C 127 : C*** LECTURE OU CALCUL DES TABLEAUX D'INDIRECTION *** 128 : C 129 : C*** On essaie de lire les tableaux d'indirection crees par PROGCS 130 : C Si la table MATRIS n'existe pas, alors LEKTAB fera appel 131 : C a PROGCS. 132 : 132 : CALL LEKTAB(MTABP,'MATRIS',MTAB2) 133 : IF(MTAB2.EQ.0)GO TO 90 134 : TYPE=' ' 134 : CALL ACMO(MTABP,'METHODE',TYPE,MTAB3) 135 : IF(MTAB3.EQ.0)GO TO 90 135 : CALL ACME(MTAB3,'KTYPI',KTYPI) 135 : CALL ACME(MTAB3,'KSTOCK',KSTO) 136 : 137 : C 138 : C*** CALCUL DES COEFFICIENTS DE LA MATRICE *** 139 : C 140 : 141 : TYPE=' ' 141 : CALL ACMO(MTABP,'MATC',TYPE,MATRAK) 142 : IF(TYPE.NE.'MATRAK')THEN 143 : WRITE(6,*)' Il n''y a pas d''entree MATC dans la table ' 144 : GO TO 90 145 : ENDIF 146 : 147 : TYPE=' ' 147 : CALL ACMO(MTABP,'DIAGV',TYPE,IZDV) 148 : IF(TYPE.NE.'CHPOINT ')THEN 149 : WRITE(6,*)' Il n''y a pas d''entree DIAGV dans la table ' 150 : GO TO 90 151 : ENDIF 152 : 153 : SEGACT MATRAK*MOD 154 : 155 : IF(KDIAG.EQ.IZDV.AND.KCAC.EQ.1)THEN 156 : GO TO 500 157 : ENDIF 158 : 159 : IF(KDPDQ.NE.0)THEN 160 : C DPDQ 161 : TYPE=' ' 161 : CALL ACMO(MTABP,'DOMAINE',TYPE,MTABD) 162 : IF(TYPE.NE.'TABLE')GO TO 90 163 : TYPE=' ' 163 : CALL ACMO(MTABD,'MELSTB',TYPE,MELSTB) 164 : IF(TYPE.NE.'MAILLAGE') GO TO 90 164 : CALL LEKTAB(MTABD,'CENTRE',MELEMC) 165 : 166 : TYPE=' ' 166 : CALL ACMO(MTABD,'MCHPOC',TYPE,MCHPOC) 167 : IF(TYPE.NE.'CHPOINT') GO TO 90 167 : CALL LICHT(MCHPOC,IZPOC,TYPC,IGEOMC) 168 : NBK=IZPOC.VPOCHA(/2) 169 : SEGACT MELEMC,IGEOMC 170 : C DPDQ 171 : ENDIF 172 : 173 : KDIAG=IZDV 173 : CALL LICHT(IZDV,MPOVAL,TYPC,IGEOM) 174 : 175 : 176 : 177 : WRITE(6,*)' *------------------------------------------------*' 178 : WRITE(6,*)' * KTYPI=2,3,4 RESOLUTION DE L''EQUATION DE *' 179 : WRITE(6,*)' * PRESSION PAR UNE METHODE ITERATIVE. *' 180 : WRITE(6,*)' * LA MATRICE EST CALCULEE ET STOCKEE UNE FOIS *' 181 : WRITE(6,*)' * POUR TOUTES. L''INVERSION UTILISE UNE *' 182 : WRITE(6,*)' * METHODE DE GRADIENT CONJUGUE PRECONDITIONNE *' 183 : WRITE(6,*)' * PAR LA DIAGONALE OU UNE FACTORISATION *' 184 : WRITE(6,*)' * INCOMPLETE. *' 185 : C 186 : C********************************************************************* 187 : C--- Initialisation des segments IZDPR IZPIM 188 : C et stockage de leurs pointeurs dans la table METHODE. 189 : C********************************************************************* 190 : C 191 : C--- NL nombre total d'elements 192 : C 192 : CALL ACME(MTAB3,'NL',NL) 193 : JG=NL 194 : SEGINI IZDPR 194 : CALL ECMO(MTAB3,'DPR','LISTREEL',IZDPR) 195 : SEGDES IZDPR 196 : IF(KTYPI.EQ.4) THEN 196 : IF(KSTO.EQ.0) CALL ACME(MTAB2,'LA',NICCG) 196 : IF(KSTO.EQ.1) CALL ACME(MTAB2,'NNZ',NNZ) 197 : IF(KSTO.EQ.1) NICCG=NL*NNZ 198 : JG=NICCG 199 : SEGINI IZICCG 199 : CALL ECMO(MTAB3,'ICCG','LISTREEL',IZICCG) 200 : SEGDES IZICCG 201 : ENDIF 202 : C 203 : SEGINI IZPIM 203 : CALL INITD (IZPIM.PROG,NL,0.D0) 203 : CALL ECMO(MTAB3,'PIM','LISTREEL',IZPIM) 204 : SEGDES IZPIM 205 : C 206 : C********************************************************************* 207 : C---- P A R A M E T R E S D E C O N V E R G E N C E 208 : C********************************************************************* 209 : C 209 : CALL ACME(MTAB3,'NITMAX',NITMAX) 209 : CALL ACMF(MTAB3,'EPSI',EPI) 209 : CALL ACME(MTAB3,'NPITE',NPITE) 209 : CALL ACME(MTAB3,'NFIMPR',NFIMPR) 210 : WRITE(6,*)' * *' 211 : WRITE(6,*)' * PARAMETRES DE CONTROLE CHOISIS : *' 212 : WRITE(6,*)' * *' 213 : WRITE(6,*)' * KTYPI : ',KTYPI 214 : WRITE(6,*)' * KSTOCK : ',KSTO 215 : WRITE(6,*)' * NITMAX : ',NITMAX 216 : WRITE(6,666) EPI 217 : WRITE(6,*)' * JTPA : ',NPITE 218 : WRITE(6,*)' * FIMPR : ',NFIMPR 219 : WRITE(6,*)' * *' 220 : 666 FORMAT(1X,' * JTEPS : ',3X,E8.2) 221 : C 222 : C********************************************************************* 223 : C---- I N I T I A L I S A T I O N D E L A P R E S S I O N 224 : C********************************************************************* 225 : C 226 : 227 : C 228 : C AU CAS ON ON FOURNIT UNE PRESSION INITIALE, ON SAUVE LE CONTENU 229 : C DANS PIM. 230 : C 231 : TYPE=' ' 231 : CALL ACMO(MTAB3,'PINITIAL',TYPE,IPRESU) 232 : C 233 : IF (IPRESU.EQ.0) THEN 234 : SEGACT IZPIM*MOD 235 : DO II=1,NL 236 : IZPIM.PROG(II)=0.D0 237 : ENDDO 238 : SEGDES IZPIM 239 : ELSE 239 : CALL LICHT(IPRESU,IPRESV,TYPC,IGEOM) 240 : SEGACT IPRESV*MOD 241 : LL1=IPRESV.VPOCHA(/2) 242 : LL2=IPRESV.VPOCHA(/1) 243 : IF (LL1.NE.1.OR.LL2.NE.NL) THEN 244 : WRITE(6,*)'PB DANS CGPR : DIM DE LA PRESSION' 245 : WRITE(6,*)'DIM ATTENDUE : ',NL 246 : WRITE(6,*)'DIM RECUE : ',LL2 247 : IF(NL.LT.LL2) THEN 248 : WRITE(6,*) 'NOMBRES DE CONTRAINTES INCOMPATIBLES' 249 : STOP 250 : ELSE 251 : WRITE(6,*) 'IL Y A DES CONTRAINTES SUPPLEMENTAIRES' 252 : SEGACT IZPIM*MOD 253 : DO II=1,LL2 254 : IZPIM.PROG(II)=IPRESV.VPOCHA(II,1) 255 : ENDDO 256 : DO II=LL2+1,NL 257 : IZPIM.PROG(II)=0.D0 258 : ENDDO 259 : SEGDES IZPIM,IPRESV 260 : ENDIF 261 : ENDIF 262 : SEGACT IZPIM*MOD 263 : DO 15 II=1,NL 264 : IZPIM.PROG(II)=IPRESV.VPOCHA(II,1) 265 : 15 CONTINUE 266 : SEGDES IZPIM,IPRESU,IPRESV 267 : ENDIF 268 : C 269 : C******************************************************************** 270 : C--- A L L O C A T I O N D E L A M A T R I C E 271 : C******************************************************************** 272 : C 273 : 274 : C--- Allocation de memoire pour la matrice A et la diagonale inverse 275 : 275 : IF(KSTO.EQ.0) CALL ACME(MTAB2,'LA',LA) 275 : IF(KSTO.EQ.1) CALL ACME(MTAB2,'NNZ',NNZ) 276 : IF(KSTO.EQ.1) LA=NL*NNZ 277 : JG=LA 278 : SEGINI IZA 278 : CALL ECMO(MTAB2,'IMAT','LISTREEL',IZA) 279 : SEGINI IZD 280 : C 281 : C******************************************************************** 282 : C--- C A L C U L D E L A M A T R I C E D E P R E S S I O N 283 : C******************************************************************** 284 : C 285 : MELEME=KLEMC 286 : C 287 : C ON COMMENCE PAR ACTIVER 288 : C 289 : 290 : MELEM1=KGEOS 290 : CALL KRIPAD(MELEM1,IZIPAD) 291 : 292 : SEGACT MELEME 293 : 294 : C 295 : C Nombre de sous-objets 296 : C 297 : NBSOUS=LISOUS(/1) 298 : IF(NBSOUS.EQ.0)NBSOUS=1 299 : 300 : IES=VPOCHA(/2) 301 : NPT=VPOCHA(/1) 302 : C 303 : 304 : C--- S T O C K A G E M O R S E 305 : IF(KSTO.EQ.0) THEN 306 : TYPE=' ' 306 : CALL ACMO(MTAB2,'IA',TYPE,IA) 307 : TYPE=' ' 307 : CALL ACMO(MTAB2,'JA',TYPE,JA) 308 : SEGACT IA*MOD,JA*MOD 309 : KKA=0 310 : C 311 : C On procede sous-objet par sous-objet 312 : C 313 : 314 : C? write(6,*)' KDPDQ=',KDPDQ 315 : IF(KDPDQ.NE.0)THEN 316 : C 317 : C OPERATEUR DPDQ 318 : C 318 : CALL KRIPAD(IGEOMC,MLENTI) 319 : SEGACT MELSTB*MOD,IGEOMC 320 : MELEMK=KGEOC 320 : CALL KRIPAD(MELEMK,IZIPAK) 321 : SEGACT MELEMK*MOD 322 : ENDIF 323 : 324 : KJ=0 325 : DO 10 JIS=1,NBSOUS 326 : IF(NBSOUS.EQ.1) IPTJ=MELEME 327 : IF(NBSOUS.NE.1) IPTJ=LISOUS(JIS) 328 : SEGACT IPTJ*MOD 329 : IPMJ=LIZAFM(JIS) 330 : JKAM=IKAM0 (JIS) 331 : SEGACT IPMJ*MOD 332 : NPJ=IPTJ.NUM(/1) 333 : NEL=IPTJ.NUM(/2) 334 : C Un sous-objet est active, on traite tous ses elements 335 : DO 1 KJJ=1,NEL 336 : KJ=KJ+1 337 : C 338 : C--- Boucle sur le nombe de voisins, les numeros de ceux-ci sont 339 : C dans JA 340 : C 341 : C On v‚rifie que la somme des termes d'une ligne de la matrice est nulle 342 : C en accumulant dans la variable STERM. 343 : C 344 : STERM=0.D0 345 : DO 2 KV=IA.LECT(KJ),IA.LECT(KJ+1)-1 346 : KKA=KKA+1 347 : KK=JA.LECT(KV) 348 : C De deux choses l'une, KK appartient au meme sous-objet 349 : C que KJ ou non. Par defaut on dit oui puis on teste. 350 : IPTK=IPTJ 351 : IPMK=IPMJ 352 : KKAM=JKAM 353 : NPK=NPJ 354 : NSOBJ=IMEM(KK) 355 : IF(NSOBJ.NE.JIS) THEN 356 : C KK et KJ sont dans deux sous-objets differents, il faut 357 : C activer celui de KK 358 : IPTK=LISOUS(NSOBJ) 359 : IPMK=LIZAFM(NSOBJ) 360 : KKAM=IKAM0 (NSOBJ) 361 : SEGACT IPTK*MOD 362 : SEGACT IPMK*MOD 363 : NPK=IPTK.NUM(/1) 364 : ENDIF 365 : KKK=KK-KKAM+1 366 : 367 : AUX=0.D0 368 : 369 : DO 31 I=1,NPJ 370 : DO 32 N=1,NPK 371 : IF(IPTJ.NUM(I,KJJ).NE.IPTK.NUM(N,KKK))GO TO 3 372 : JD=IZIPAD.LECT(IPTK.NUM(N,KKK)) 373 : * IPR=IC*(JD-1)+1 374 : AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1) 375 : & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2) 376 : IF(IES.EQ.2)GO TO 3 377 : AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3) 378 : 3 CONTINUE 379 : 32 CONTINUE 380 : 31 CONTINUE 381 : IZA.PROG(KKA)=AUX 382 : IF(KK.EQ.KJ) IZD.PROG(KJ)=AUX 383 : 384 : C? write(6,*)' KDPDQ 2 = ',KDPDQ 385 : IF(KDPDQ.NE.0)THEN 386 : C 387 : C OPERATEUR DPDQ 388 : C 389 : C PV nuna est dans le segment izl. izl n'est pas initialisé 389 : if (izl.eq.0) call erreur(5) 390 : KJQ=NUNA(KJ) 391 : LPJ=LECT(MELEMK.NUM(1,KJQ)) 392 : IF(LPJ.NE.0)THEN 393 : P1=BETA*IZPOC.VPOCHA(LPJ,1) 394 : 395 : C?! D(KJ)=D(KJ)+ABS(P1) 396 : IZD.PROG(KJ)=IZD.PROG(KJ)+ABS(P1) 397 : 398 : DO 3819 JC=2,NBK 399 : NPC=MELSTB.NUM(JC,LPJ) 400 : P=BETA*IZPOC.VPOCHA(LPJ,JC) 401 : KPI=IZIPAK.LECT(NPC) 402 : KPJN=NUAN(KPI) 403 : IF(KPJN.GT.KJ)GO TO 3819 404 : IF(KPJN.GT.KJ)GO TO 3819 405 : 406 : C# MC : IDEBLK non initialisé + IDECI et LLI uniquement dans KJAA 407 : C KJAA non utilisé => on vire !!! 408 : C IDECI=IDEBLK(KJ-KJD+1) 409 : C LLI=IDEBLK(KJ-KJD+2)-IDEBLK(KJ-KJD+1) 410 : C KJAA=LLI-KJ+KPJN+IDECI 411 : 412 : C?! A(KJAA)=A(KJAA)+P 413 : IZA.PROG(KKA)=IZA.PROG(KKA)+P 414 : 415 : 3819 CONTINUE 416 : ENDIF 417 : ENDIF 418 : C FIN DPDQ 419 : 420 : C 421 : STERM=STERM+AUX 422 : C 423 : C Il faut desactiver IPMK et IPTK sauf bien sur si on n'a 424 : C qu'un seul sous-objet 425 : C 426 : IF(NSOBJ.NE.JIS) THEN 427 : SEGDES IPMK 428 : SEGDES IPTK 429 : ENDIF 430 : 2 CONTINUE 431 : C 432 : C Controle sur la somme des termes d'une ligne 433 : C 434 : IF(SOMTER.AND.ABS(STERM).GT.1.D-10) THEN 435 : WRITE(6,*) 'L''ELEMENT KJ = ',KJ,' DOIT ETRE SUR UNE SORTIE' 436 : ENDIF 437 : 1 CONTINUE 438 : SEGDES IPTJ 439 : SEGDES IPMJ 440 : 10 CONTINUE 441 : C 442 : C Sortie 443 : C dimension ts(9) 444 : C do 1001 kj=1,nelp 445 : C write(61,60) (itab(ka),ka=iat(kj),iat(kj+1)-1) 446 : C1001 continue 447 : C do 1000 kj=1,nelp 447 : C call initd(ts,9,0.d0) 448 : C do 1003 ka=iat(kj),iat(kj+1)-1 449 : C i=ka-iat(kj)+1 450 : C ts(i)=a(ka) 451 : C1003 continue 452 : C write(61,61) (a(ka),ka=iat(kj),iat(kj+1)-1) 453 : C write(61,61) (ts(i),i=1,9) 454 : C1000 continue 455 : 60 format(9(1x,i4)) 456 : 61 format(9(1x,f7.4)) 457 : C 458 : SEGDES IA,JA,IZA 459 : SEGDES MELEME 460 : 461 : IF(KDPDQ.NE.0)THEN 462 : SEGSUP IZIPAK,MLENTI 463 : ENDIF 464 : 465 : 466 : C--- S T O C K A G E C O M P R E S S E 467 : 468 : ELSEIF(KSTO.EQ.1) THEN 469 : TYPE=' ' 469 : CALL ACMO(MTAB2,'KA',TYPE,KA) 470 : SEGACT KA*MOD 471 : KJ=0 472 : DO 20 JIS=1,NBSOUS 473 : IF(NBSOUS.EQ.1) IPTJ=MELEME 474 : IF(NBSOUS.NE.1) IPTJ=LISOUS(JIS) 475 : SEGACT IPTJ*MOD 476 : IPMJ=LIZAFM(JIS) 477 : JKAM=IKAM0 (JIS) 478 : SEGACT IPMJ*MOD 479 : NPJ=IPTJ.NUM(/1) 480 : NEL=IPTJ.NUM(/2) 481 : C Un sous-objet est active, on traite tous ses elements 482 : DO 101 KJJ=1,NEL 483 : KJ=KJ+1 484 : C 485 : C--- Boucle sur le nombre de voisins, la liste se termine lorsque 486 : C le numero du voisin est egal au numero de l'element. 487 : C KA et A sont ranges comme si ils etaient declares KA(NL,NNZ) 488 : C et A(NL,NNZ). 489 : C 490 : STERM=0.D0 491 : DO 102 KV=1,NNZ 492 : KKA=KJ+(KV-1)*NL 493 : KK=KA.LECT(KKA) 494 : IF(KV.GT.1.AND.KK.EQ.KJ) THEN 495 : C--- On est arrive a la fin de la liste des voisins 496 : IZA.PROG(KKA)=0.D0 497 : GOTO 102 498 : ENDIF 499 : C De deux choses l'une, KK appartient au meme sous-objet 500 : C que KJ ou non. Par defaut on dit oui puis on teste. 501 : IPTK=IPTJ 502 : IPMK=IPMJ 503 : KKAM=JKAM 504 : NPK=NPJ 505 : NSOBJ=IMEM(KK) 506 : IF(NSOBJ.NE.JIS) THEN 507 : C KK et KJ sont dans deux sous-objets differents, il faut 508 : C activer celui de KK 509 : IPTK=LISOUS(NSOBJ) 510 : IPMK=LIZAFM(NSOBJ) 511 : KKAM=IKAM0 (NSOBJ) 512 : SEGACT IPTK*MOD 513 : SEGACT IPMK*MOD 514 : NPK=IPTK.NUM(/1) 515 : ENDIF 516 : KKK=KK-KKAM+1 517 : 518 : AUX=0.D0 519 : DO 131 I=1,NPJ 520 : DO 132 N=1,NPK 521 : IF(IPTJ.NUM(I,KJJ).NE.IPTK.NUM(N,KKK))GO TO 103 522 : JD=IZIPAD.LECT(IPTK.NUM(N,KKK)) 523 : * IPR=IC*(JD-1)+1 524 : AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1) 525 : & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2) 526 : IF(IES.EQ.2)GO TO 103 527 : AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3) 528 : 103 CONTINUE 529 : 132 CONTINUE 530 : 131 CONTINUE 531 : IZA.PROG(KKA)=AUX 532 : IF(KV.EQ.1) IZD.PROG(KJ)=AUX 533 : C 534 : C Il faut desactiver IPMK et IPTK sauf bien sur si on n'a 535 : C qu'un seul sous-objet 536 : C 537 : STERM=STERM+AUX 538 : IF(NSOBJ.NE.JIS) THEN 539 : SEGDES IPMK 540 : SEGDES IPTK 541 : ENDIF 542 : 102 CONTINUE 543 : IF(SOMTER.AND.ABS(STERM).GT.1.D-10) THEN 544 : WRITE(6,*) 'L''ELEMENT KJ = ',KJ,' DOIT ETRE SUR UNE SORTIE' 545 : ENDIF 546 : 101 CONTINUE 547 : SEGDES IPTJ 548 : SEGDES IPMJ 549 : 20 CONTINUE 550 : C 551 : C Sortie 552 : C do 2001 kj=1,nelp 553 : C write(60,60) (itab(kj+(kv-1)*nelp),kv=1,nnz) 554 : C2001 continue 555 : C do 2000 kj=1,nelp 556 : C write(60,61) ( a (kj+(kv-1)*nelp),kv=1,nnz) 557 : C2000 continue 558 : SEGDES KA,IZA 559 : SEGDES MELEME 560 : ELSE 561 : WRITE(6,*) 'CGPR Mode de stockage inconnu |' 562 : ENDIF 563 : C 564 : C*** RAJOUTER LA DIAGONALE DE LA PRESSION *** 565 : C 566 : TYPE=' ' 566 : CALL ACMO(MTABP,'DIAGP',TYPE,IZDP) 567 : IF(TYPE.EQ.'CHPOINT ')THEN 567 : CALL LICHT(IZDP,IZDDP,TYPE,IGEOMI) 568 : IGEOM=KGEOC 569 : 570 : SEGACT IZD*MOD 571 : SEGACT IGEOM,IGEOMI 572 : 573 : NPT=IZD.PROG(/1) 574 : NPTI=IZDDP.VPOCHA(/1) 575 : 576 : IF(KSTO.EQ.0) THEN 577 : SEGACT IA*MOD,JA*MOD,IZA*MOD 578 : ELSE 579 : SEGACT KA*MOD,IZA*MOD 580 : ENDIF 581 : DO 51 I=1,NPTI 582 : XVALI=IZDDP.VPOCHA(I,1) 583 : IPOI=IGEOMI.NUM(1,I) 584 : DO 52 J=1,NPT 585 : IF(IGEOM.NUM(1,J).EQ.IPOI)THEN 586 : IZD.PROG(J)=IZD.PROG(J)+XVALI 587 : C 588 : C*** Il faut aussi mettre a jour la matrice 589 : C 590 : IF(KSTO.EQ.0) THEN 591 : IZA.PROG(IA.LECT(J))=IZA.PROG(IA.LECT(J))+XVALI 592 : ELSE 593 : IZA.PROG(J)=IZA.PROG(J)+XVALI 594 : ENDIF 595 : 596 : GO TO 53 597 : ENDIF 598 : 52 CONTINUE 599 : 53 CONTINUE 600 : 51 CONTINUE 601 : IF(KSTO.EQ.0) THEN 602 : SEGDES IA,JA,IZA 603 : ELSE 604 : SEGDES KA,IZA 605 : ENDIF 606 : SEGDES IZDP,IZDDP,IGEOM,IGEOMI 607 : ENDIF 608 : C 609 : C******************************************************************** 610 : C--- P R E C O N D I T I O N N E M E N T 611 : C******************************************************************** 612 : C 613 : 614 : C write(6,*)' P R E C O N' 615 : IF(KTYPI.EQ.2) THEN 616 : SEGACT IZDPR*MOD 616 : CALL INITD(IZDPR.PROG,NL,1.D0) 616 : CALL ECMO(MTAB3,'IPRC','LISTREEL',IZDPR) 617 : SEGDES IZDPR 618 : ELSEIF(KTYPI.EQ.3) THEN 619 : SEGACT IZDPR*MOD 620 : DO 23 K=1,NL 621 : IZDPR.PROG(K)=1./SQRT(IZD.PROG(K)) 622 : 23 CONTINUE 622 : CALL ECMO(MTAB3,'IPRC','LISTREEL',IZDPR) 623 : SEGDES IZDPR 624 : ELSEIF(KTYPI.EQ.4) THEN 625 : IF(KSTO.EQ.0) THEN 626 : SEGACT IZICCG*MOD,IZA*MOD,IA*MOD,JA*MOD 626 : CALL ICGFAM(IZICCG.PROG,IZA.PROG,IA.LECT,JA.LECT,NL) 627 : SEGDES IZICCG,IZA,IA,JA 628 : ENDIF 629 : IF(KSTO.EQ.1) THEN 630 : SEGACT IZICCG*MOD,IZA*MOD,KA*MOD 630 : CALL ICGFAC(IZICCG.PROG,IZA.PROG,KA.LECT,NL,NNZ) 631 : SEGDES IZICCG,IZA,KA 632 : ENDIF 632 : CALL ECMO(MTAB3,'IPRC','LISTREEL',IZICCG) 633 : ELSE 634 : WRITE(6,*) 'CGPR Option inconnue' 635 : END IF 636 : SEGSUP IZD 637 : SEGSUP IZIPAD 638 : SEGDES MPOVAL,IZDV 639 : KCAC=1 640 : SEGDES MATRAK 641 : SEGDES MELEME 642 : C 643 : C*** RESOLUTION *** 644 : C 645 : 500 CONTINUE 645 : CALL LICHT(MCHB,IZB,TYPC,IGEOM) 645 : CALL LICHT(MCHR,IZR,TYPC,IGEOM) 646 : C write(6,*)' AVT KAL2P ' 646 : CALL KAL2P(MTABP,MCHB) 647 : C write(6,*)' AVT RESCGS' 647 : CALL RESCGS(IZB,MTABP,IZR) 648 : C write(6,*)' APR RESCGS' 649 : segact melstb 650 : RETURN 651 : 90 CONTINUE 652 : WRITE(6,*)' Arret anormal de KRESS' 653 : RETURN 654 : END 655 : 656 : 657 : 658 : 659 : 660 : 661 : 662 :
© Cast3M 2003 - All rights reserved.
Disclaimer