1 : C KMDMT SOURCE CHAT 06/06/01 21:16:56 5450 2 : SUBROUTINE KMDMT(MTABP,MCHB,MCHR,IMPR,CBETA,KDPDQ,KPIMP,PIMP,NIMP) 3 : C *********************************************************************** 4 : C -1 T 5 : C CE SOUS PROGRAMME CALCULE LA MATRICE C D C CONNAISSANT C 6 : C 7 : C POINTEURS : 8 : C 9 : C IZL CONTIENT LE PROFIL ET LA RENUMEROTATION (CF PROFCH) 10 : C MELEME OBJET MAILLAGE SUR LEQUEL EST FAIT LE CALCUL 11 : C MATRAK MATRICES ELEMENTAIRES DE LA DIVERGENCE (ALIAS"C") 12 : C 13 : C IZDV DIAGONALE D 14 : C IZIPAD CORRESPONDANCE NUMER. GLOBALE --> NUMER. LOCALE 15 : C (DOMAINE SUR LEQUEL PORTE AP ET NON LA PRESSION) 16 : C 17 : C EN SORTIE : 18 : C 19 : C LA DEMI-MATRICE INFERIEUR EST RANGEE PAR LIGNE 20 : C (UN SEGMENT IZA PAR LIGNE) 21 : C LA DIAGONALE DE LA MATRICE EST RANGEE A PART DANS IZD 22 : C KZA CONTIENT TOUT LES POINTEURS IZA 23 : C IZD=KZA(1) 24 : C 25 : C *********************************************************************** 26 : IMPLICIT INTEGER(I-N) 27 : IMPLICIT REAL*8 (A-H,O-Z) 28 : 29 : -INC CCOPTIO 30 : -INC SMELEME 31 : POINTEUR IPTJ.MELEME,IPTK.MELEME,MELEM1.MELEME 32 : POINTEUR IGEOM.MELEME,IGEOMI.MELEME,IGEOMC.MELEME 33 : POINTEUR MELEMK.MELEME,MELEMC.MELEME 34 : POINTEUR MELSTB.MELEME 35 : C? POINTEUR ELTFA.MELEME,FACEL.MELEME,MELEMC.MELEME 36 : C? POINTEUR MELEMF.MELEME,MELEMK.MELEME,MELELI.MELEME 37 : -INC SMLENTI 38 : POINTEUR IZIPAD.MLENTI,MLENTP.MLENTI,IZIPAK.MLENTI 39 : -INC SMCOORD 40 : -INC SMCHPOI 41 : POINTEUR IZDV.MCHPOI 42 : POINTEUR IZDP.MCHPOI,IZDDP.MPOVAL 43 : POINTEUR MCHB.MCHPOI,MCHR.MCHPOI 44 : C? POINTEUR IZB.MPOVAL,IZR.MPOVAL,IZPOF.MPOVAL,IZPOC.MPOVAL 45 : POINTEUR IZB.MPOVAL,IZR.MPOVAL,IZPOC.MPOVAL 46 : C? POINTEUR IZP.MPOVAL,IZF.MPOVAL 47 : POINTEUR IZP.MPOVAL 48 : C? POINTEUR IZFAC.MPOVAL,IZVOL.MPOVAL 49 : 50 : C-INC SMMATRAKANC 51 : -INC SMMATRK1 52 : C************************************************************************* 53 : C 54 : C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees 55 : C 56 : 57 : * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange 58 : * (points CENTRE ) pour chaque operateur de contrainte 59 : * KGEOC SPG pour la totalite des points CENTRE. 60 : * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse) 61 : * KLEMC Connectivites de l'ensemble des contraintes 62 : * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones 63 : 64 : SEGMENT MATRAK 65 : INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP) 66 : INTEGER LIZAFM(NBSOUS) 67 : INTEGER IKAM0 (NBSOUS) 68 : INTEGER IMEM (NBELC) 69 : INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC 70 : ENDSEGMENT 71 : 72 : SEGMENT IZAFM 73 : REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX) 74 : ENDSEGMENT 75 : 76 : POINTEUR IPMJ.IZAFM,IPMK.IZAFM 77 : 78 : C******************************************************************* 79 : SEGMENT/IZA/(A(LA)*D) 80 : SEGMENT/IZD/(D(NELP)*D) 81 : -INC SMTABLE 82 : POINTEUR MTABP.MTABLE 83 : C 84 : CHARACTER*8 TYPE,TYPC 85 : * SAVE BETA 86 : C*** 87 : C KDPDQ=0 88 : C write(6,*)' Entree dans kmdmt KDPDQ=',kdpdq 89 : 90 : TYPE=' ' 90 : CALL ACMO(MTABP,'MATC',TYPE,MATRAK) 91 : C WRITE(6,*)' KMDMT : MATC=',MATRAK 92 : IF(TYPE.NE.'MATRAK')THEN 93 : WRITE(6,*)' IL N''Y A PAS D''ENTREE MATC DANS LA TABLE ' 94 : GO TO 90 95 : ENDIF 96 : 97 : SEGACT MATRAK 98 : * initialisations pour plaire a l'optimiseur IBM 99 : MELEMK=MATRAK 100 : IZPOC=MATRAK 101 : IZIPAK=MATRAK 102 : MELSTB=MATRAK 103 : IZA=MATRAK 104 : IZL=MATRAK 105 : MLENTI=MATRAK 106 : IZD=MATRAK 107 : IDBLK=MATRAK 108 : IPMJ=MATRAK 109 : NBELC=IMEM(/1) 110 : NIMP=NBELC/2 111 : IF(KPIMP.EQ.2)NIMP=NBELC 112 : C write(6,*)' CAS KPIMP=',kpimp,' PIMP=',pimp 113 : 113 : CALL PROFCH(MATRAK,IMPR,KPIMP,NIMP) 114 : 115 : TYPE=' ' 115 : CALL ACMO(MTABP,'DIAGV',TYPE,IZDV) 116 : IF(TYPE.NE.'CHPOINT ')THEN 117 : WRITE(6,*)' IL N''Y A PAS D''ENTREE DIAGV DANS LA TABLE ' 118 : GO TO 90 119 : ENDIF 120 : 121 : SEGACT MATRAK*MOD 122 : IZL=KIZCL 123 : C write(6,*)' KIZCL=',KIZCL 124 : C WRITE(6,*)' KDIAG=',KDIAG,' IZDV=',IZDV,' KCAC=',KCAC 125 : IF(KDIAG.EQ.IZDV.AND.KCAC.EQ.1)THEN 126 : GO TO 500 127 : ENDIF 128 : 129 : 130 : IF(KDPDQ.NE.0)THEN 131 : C DPDQ 132 : TYPE=' ' 132 : CALL ACMO(MTABP,'DOMAINE',TYPE,MTABD) 133 : IF(TYPE.NE.'TABLE')GO TO 90 134 : TYPE=' ' 134 : CALL ACMO(MTABD,'MELSTB',TYPE,MELSTB) 135 : IF(TYPE.NE.'MAILLAGE') GO TO 90 135 : CALL LEKTAB(MTABD,'CENTRE',MELEMC) 136 : 137 : TYPE=' ' 137 : CALL ACMO(MTABD,'MCHPOC',TYPE,MCHPOC) 138 : IF(TYPE.NE.'CHPOINT') GO TO 90 138 : CALL LICHT(MCHPOC,IZPOC,TYPC,IGEOMC) 139 : NBK=IZPOC.VPOCHA(/2) 140 : SEGACT MELEMC,IGEOMC 141 : C DPDQ 142 : ENDIF 143 : 144 : KDIAG=IZDV 144 : CALL LICHT(IZDV,MPOVAL,TYPC,IGEOM) 145 : 146 : C WRITE(6,*)' DIAGONALE : ' 147 : C LDG=VPOCHA(/1) 148 : C WRITE(6,*)' SUR X ' 149 : C WRITE(6,1002)(VPOCHA(II,1),II=1,LDG) 150 : C WRITE(6,*)' SUR Y ' 151 : C WRITE(6,1002)(VPOCHA(II,2),II=1,LDG) 152 : 153 : C 154 : C ON COMMENCE PAR ACTIVER 155 : C 156 : SEGACT IZL 157 : 158 : C. NJAN=NUAN(/1) 159 : C. write(6,*)' NUAN=' 160 : C. write(6,1001)(NUAN(II),II=1,NJAN) 161 : C. write(6,*)' NUNA=' 162 : C. write(6,1001)(NUNA(II),II=1,NJAN) 163 : C write(6,*)' IMEL=' 164 : C write(6,1001)(IMEL(II),II=1,NJAN) 165 : C write(6,*)' IMJ=' 166 : C write(6,1001)(IMJ(II),II=1,NJAN) 167 : 168 : MELEM1=KGEOS 168 : CALL KRIPAD(MELEM1,IZIPAD) 169 : 170 : MELEME=KLEMC 171 : SEGACT MELEME 172 : 173 : IKAS=0 174 : IF(IKAS.EQ.1)THEN 175 : NIZ=LIZAFM(/1) 176 : DO 1940 L=1,NIZ 177 : IZAFM=LIZAFM(L) 178 : WRITE(6,*)' MATRICE AM BLOC N° ',L 179 : SEGACT IZAFM 180 : NP=AM(/2) 181 : IES=AM(/3) 182 : NNELP=AM(/1) 183 : WRITE(6,*)'NP,IES,NNELP=',NP,IES,NNELP 184 : DO 1941 K=1,NNELP 185 : WRITE(6,*)'K=',K 186 : WRITE(6,1008)((AM(K,I,J),I=1,NP),J=1,IES) 187 : 1941 CONTINUE 188 : 1940 CONTINUE 188 : C CALL ARRET(0) 189 : 1008 FORMAT(8(1X,1PE11.4)) 190 : ENDIF 191 : 192 : IES=VPOCHA(/2) 193 : NPT=VPOCHA(/1) 194 : NELP=B(/1) 195 : C 196 : IDMAT=KZA1 197 : SEGACT IDMAT*MOD 198 : NBLK=IDESCR(/1) 199 : WRITE(6,*) '* NOMBRE DE BLOCS DE LA MATRICE = ',NBLK 200 : SEGINI IZD 201 : IDIAG=IZD 202 : C 203 : 204 : DO 71 NS=1,LIZAFM(/1) 205 : IZAFM=LIZAFM(NS) 206 : SEGACT IZAFM 207 : C 208 : C NP1=AM(/1) 209 : C NP2=AM(/3) 210 : C WRITE(6,*)' NP=',NP1,' NELB=',NP2 211 : C DO 1834 I=1,NP1 212 : C WRITE(6,*)' I1=',I 213 : C WRITE(6,1002)(AM(I,1,KK),KK=1,NP2) 214 : C WRITE(6,*)' I2=',I 215 : C WRITE(6,1002)(AM(I,2,KK),KK=1,NP2) 216 : C1834 CONTINUE 217 : 218 : 71 CONTINUE 219 : 220 : NBSOUS=LISOUS(/1) 221 : IF(NBSOUS.EQ.0)NBSOUS=1 222 : 223 : DO 10 JIS=1,NBSOUS 224 : IF(NBSOUS.EQ.1)IPTJ=MELEME 225 : IF(NBSOUS.NE.1)IPTJ=LISOUS(JIS) 226 : SEGACT IPTJ 227 : 10 CONTINUE 228 : 229 : IF(KDPDQ.NE.0)THEN 230 : C 231 : C OPERATEUR DPDQ 232 : C 232 : CALL KRIPAD(IGEOMC,MLENTI) 233 : SEGACT MELSTB,IGEOMC 234 : MELEMK=KGEOC 234 : CALL KRIPAD(MELEMK,IZIPAK) 235 : SEGACT MELEMK 236 : ENDIF 237 : 238 : C Calcul de la diagonale 239 : 240 : AD=0.D0 241 : ADB=1.D0 242 : NBELP=NUNA(/1) 243 : DO 7710 KJ=1,NBELP 244 : KJJ=NUNA(KJ) 245 : 246 : IPMJ=LIZAFM(IMEM(NUNA(KJ))) 247 : JKAM=IKAM0 (IMEM(NUNA(KJ))) 248 : KJJ=KJJ-JKAM+1 249 : J=IMJ(KJ) 250 : IPTJ=IMEL(KJ) 251 : NPJ=IPTJ.NUM(/1) 252 : K1=KJ-KZA(KJ)+1 253 : 254 : KK=KJ 255 : KKK=NUNA(KK) 256 : 257 : IPMK=LIZAFM(IMEM(NUNA(KK))) 258 : KKAM=IKAM0 (IMEM(NUNA(KK))) 259 : KKK=KKK-KKAM+1 260 : K=IMJ(KK) 261 : IPTK=IMEL(KK) 262 : NPK=IPTK.NUM(/1) 263 : IF(KK.NE.KJ) KA=KA+1 264 : AUX=0.D0 265 : DO 7731 I=1,NPJ 266 : DO 7732 N=1,NPK 267 : IF(IPTJ.NUM(I,J).NE.IPTK.NUM(N,K))GO TO 7732 268 : JD=IZIPAD.LECT(IPTK.NUM(N,K)) 269 : 270 : AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1) 271 : & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2) 272 : IF(IES.EQ.3) 273 : &AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3) 274 : 275 : 7732 CONTINUE 276 : 7731 CONTINUE 277 : C 278 : AD=AD+AUX 279 : D(KJ)=D(KJ)+AUX 280 : 281 : 7710 CONTINUE 282 : AD=AD/FLOAT(NBELP) 283 : 284 : 285 : IF(KDPDQ.NE.0)THEN 286 : C 287 : C OPERATEUR DPDQ 288 : C 289 : NBELP=NUNA(/1) 290 : ADB=0.D0 291 : DO 7711 KJ=1,NBELP 292 : KC=NUNA(KJ) 293 : LPC=LECT(MELEMK.NUM(1,KC)) 294 : IF(LPC.NE.0)THEN 295 : ADB=ADB+ABS(IZPOC.VPOCHA(LPC,1)) 296 : ENDIF 297 : 7711 CONTINUE 298 : ADB=ADB/FLOAT(NBELP) 299 : 300 : BETA0=AD/ADB*5.E-4 301 : BETA=CBETA*BETA0 302 : C WRITE(6,7271)AD,ADB,BETA0,CBETA 303 : C7271 FORMAT(1X,' AD=',1PE11.4,' ADB=',1PE11.4,' BETAopt=',1PE11.4, 304 : WRITE(6,7271)BETA0,CBETA 305 : 7271 FORMAT(1X,' BETAopt=',1PE11.4,' CBETA=',1PE11.4) 306 : 307 : KJ=1 308 : KC=NUNA(KJ) 309 : LPC=LECT(MELEMK.NUM(1,KC)) 310 : IF(LPC.NE.0)THEN 311 : P1=BETA*IZPOC.VPOCHA(LPC,1) 312 : D(KJ)=D(KJ)+ABS(P1) 313 : ENDIF 314 : 315 : ENDIF 316 : C FIN DPDQ 317 : 318 : C 319 : C CALCUL DES BLOCS 320 : C 321 : DO 100 IBLK=1,NBLK 322 : C NUMÉROS DES LIGNES DE DÉBUT ET DE FIN DE BLOC KJD,KJF 323 : KJD=NLDBLK(IBLK) 324 : KJF=NLDBLK(IBLK+1)-1 325 : C ALLOCATION DE LA MÉMOIRE POUR LE BLOC 326 : IDBLK=IDESCR(IBLK) 327 : SEGACT IDBLK*MOD 328 : LA=ILON 329 : WRITE(6,*) '* BLOC NUMÉRO : ',IBLK,' TAILLE = ',ILON,' MOTS' 330 : SEGINI IZA 331 : IMAT=IZA 332 : KA=0 333 : DO 1 KJ=KJD,KJF 334 : KJJ=NUNA(KJ) 335 : 336 : IPMJ=LIZAFM(IMEM(NUNA(KJ))) 337 : JKAM=IKAM0 (IMEM(NUNA(KJ))) 338 : KJJ=KJJ-JKAM+1 339 : J=IMJ(KJ) 340 : IPTJ=IMEL(KJ) 341 : NPJ=IPTJ.NUM(/1) 342 : K1=KJ-KZA(KJ)+1 343 : 344 : DO 2 KK=K1,KJ 345 : KKK=NUNA(KK) 346 : 347 : IPMK=LIZAFM(IMEM(NUNA(KK))) 348 : KKAM=IKAM0 (IMEM(NUNA(KK))) 349 : KKK=KKK-KKAM+1 350 : K=IMJ(KK) 351 : IPTK=IMEL(KK) 352 : NPK=IPTK.NUM(/1) 353 : IF(KK.NE.KJ) KA=KA+1 354 : AUX=0.D0 355 : DO 31 I=1,NPJ 356 : DO 32 N=1,NPK 357 : IF(IPTJ.NUM(I,J).NE.IPTK.NUM(N,K))GO TO 32 358 : JD=IZIPAD.LECT(IPTK.NUM(N,K)) 359 : 360 : AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1) 361 : & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2) 362 : IF(IES.EQ.3) 363 : &AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3) 364 : 365 : 32 CONTINUE 366 : 31 CONTINUE 367 : C 368 : IF(KK.EQ.KJ)GO TO 21 369 : A(KA)=A(KA)+AUX 370 : C 371 : 21 CONTINUE 372 : 2 CONTINUE 373 : 374 : C?. ad=ad+aux 375 : C?. D(KJ)=D(KJ)+AUX 376 : 377 : IF(KDPDQ.NE.0)THEN 378 : C OPERATEUR DPDQ 379 : C 380 : KJJ=NUNA(KJ) 381 : LPJ=LECT(MELEMK.NUM(1,KJJ)) 382 : IF(LPJ.NE.0)THEN 383 : P1=BETA*IZPOC.VPOCHA(LPJ,1) 384 : D(KJ)=D(KJ)+ABS(P1) 385 : 386 : DO 3819 JC=2,NBK 387 : NPC=MELSTB.NUM(JC,LPJ) 388 : P=BETA*IZPOC.VPOCHA(LPJ,JC) 389 : KPI=IZIPAK.LECT(NPC) 390 : KPJN=NUAN(KPI) 391 : IF(KPJN.GT.KJ)GO TO 3819 392 : IDECI=IDEBLK(KJ-KJD+1) 393 : LLI=IDEBLK(KJ-KJD+2)-IDEBLK(KJ-KJD+1) 394 : KJAA=LLI-KJ+KPJN+IDECI 395 : 396 : A(KJAA)=A(KJAA)+P 397 : 398 : 3819 CONTINUE 399 : ENDIF 400 : ENDIF 401 : C FIN DPDQ 402 : 403 : 1 CONTINUE 404 : SEGDES IZA 405 : SEGDES IDBLK 406 : 100 CONTINUE 407 : SEGDES IDMAT 408 : 409 : 410 : IF(KDPDQ.NE.0)THEN 411 : SEGSUP IZIPAK,MLENTI 412 : ENDIF 413 : 414 : 415 : IF(KPIMP.EQ.1)THEN 416 : D(NIMP)=1.D30 417 : ELSEIF(KPIMP.EQ.2)THEN 418 : SEGACT IZA 419 : NELP=D(/1) 420 : LGLA=A(/1) 421 : LGLI=KZA(NIMP) 421 : CALL INITD(A(LGLA-LGLI+1),LGLI,0.D0) 422 : D(NELP)=1.D-10 422 : CALL LEKTAB(MTABD,'XXVOLUM',MCHPOI) 422 : CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) 423 : VOLT=0.D0 424 : DO 482 I=1,NELP 425 : VOLT=VOLT+VPOCHA(I,1) 426 : 482 CONTINUE 427 : DO 483 I=1,NELP 428 : XVOL=VPOCHA(NUNA(I),1) 429 : A(LGLA-LGLI+I)=-XVOL/VOLT 430 : 483 CONTINUE 431 : ENDIF 432 : 433 : SEGSUP IZIPAD 434 : 435 : DO 20 JIS=1,NBSOUS 436 : IF(NBSOUS.EQ.1)IPTJ=MELEME 437 : IF(NBSOUS.NE.1)IPTJ=LISOUS(JIS) 438 : C write(6,*)' On desactive JIS,IPTJ=',JIS,IPTJ 439 : SEGDES IPTJ 440 : IPMJ=LIZAFM(JIS) 441 : JKAM=IKAM0 (JIS) 442 : SEGDES IPMJ 443 : 20 CONTINUE 444 : 445 : SEGDES MPOVAL,IZDV 446 : KCAC=1 447 : 448 : TYPE=' ' 448 : CALL ACMO(MTABP,'DIAGP',TYPE,IZDP) 449 : IF(TYPE.EQ.'CHPOINT ')THEN 449 : CALL LICHT(IZDP,IZDDP,TYPE,IGEOMI) 450 : IGEOM=KGEOC 451 : 452 : SEGACT IZD 453 : SEGACT IGEOM,IGEOMI 454 : 455 : NPT=D(/1) 456 : NPTI=IZDDP.VPOCHA(/1) 457 : 458 : DO 51 I=1,NPTI 459 : XVALI=IZDDP.VPOCHA(I,1) 460 : IPOI=IGEOMI.NUM(1,I) 461 : DO 52 J=1,NPT 462 : IF(IGEOM.NUM(1,J).EQ.IPOI)THEN 463 : JJ=NUAN(J) 464 : D(JJ)=D(JJ)+XVALI 465 : GO TO 53 466 : ENDIF 467 : 52 CONTINUE 468 : 53 CONTINUE 469 : 51 CONTINUE 470 : SEGDES IZDP,IZDDP 471 : SEGDES IGEOM,IGEOMI 472 : ENDIF 473 : 474 : DO 72 NS=1,LIZAFM(/1) 475 : IZAFM=LIZAFM(NS) 476 : SEGDES IZAFM 477 : 72 CONTINUE 478 : 479 : SEGDES IZD 480 : SEGDES MATRAK 481 : SEGDES MELEME 482 : C WRITE(6,*)' APPEL A TRIAKS' 483 : 483 : CALL TRIAKS(IZL) 484 : C WRITE(6,*)' RETOUR DE TRIAKS' 485 : 500 CONTINUE 485 : CALL LICHT(MCHB,IZB,TYPC,IGEOM) 485 : CALL LICHT(MCHR,IZR,TYPC,IGEOM) 486 : 486 : CALL KAL2P(MTABP,MCHB) 487 : 488 : C write(6,*)' APPEL A RESOCK' 489 : IF(KPIMP.NE.0)THEN 490 : SEGACT IZL 490 : CALL LICHT(MCHB,IZB,TYPC,IGEOM) 491 : IF(KPIMP.EQ.1)IZB.VPOCHA(NUNA(NIMP),1)=PIMP*1.D30 492 : IF(KPIMP.EQ.2)IZB.VPOCHA(NUNA(NIMP),1)=-PIMP 493 : ENDIF 494 : 494 : CALL RESOCK(IZB,IZL,IZR) 495 : C write(6,*)' Retour de resock ' 496 : 497 : SEGDES MCHB,IZB 498 : SEGDES MCHR,IZR 499 : RETURN 500 : 90 CONTINUE 501 : WRITE(6,*)' ARRET ANORMAL DE KMDMT' 502 : RETURN 503 : 1002 FORMAT(10(1X,1PE11.4)) 504 : 1001 FORMAT(20(1X,I5)) 505 : END 506 : 507 : 508 : 509 : 510 : 511 : 512 : 513 :
© Cast3M 2003 - All rights reserved.
Disclaimer