1 : C ASSEM2 SOURCE PV 11/05/03 21:15:00 6959 2 : SUBROUTINE ASSEM2(ITRAV1,ITOPO1,INUIN1,IMINI1,MMMTRI,IPO1,INCTR1 3 : $ ,IITOP1) 4 : C 5 : C **** SUBROUTINE POUR FAIRE L'ASSEMBLAGE DE MATRICES SYMETRIQUES 6 : C EN VUE D'UN TRAITEMENT PAR METHODE DE KROUT. 7 : C 8 : C EN ENTREE: 9 : C **** ITRAV1 : POINTEUR OBJET MRIGIDITE 10 : C **** ITOPO1 : POINTEUR SEGMENT DE TRAVAIL ITOPO ( VOIR ASSEM1) 11 : C **** IITOP1 : POINTEUR SEGMENT DE TRAVAIL IITOP ( VOIR ASSEM1) 12 : C **** INUIN1 : POINTEUR SEGMENT DE TRAVAIL INUINV(VOIR ASSEM1) 13 : C **** IMINI1 : POINTEUR SEGMENT DE TRAVAIL IMINI (VOIR ASSEM1) 14 : C **** IPO1 : POINTEUR SEGMENT DE TRAVAIL IPOS (VOIR ASSEM1) 15 : C **** MMMTRI : POINTEUR OBJET MATRICE TRIANGULARISEE (NON MODIFIE) 16 : C (VOIR SMMATRI) 17 : IMPLICIT REAL*8 (A-H,O-Z) 18 : IMPLICIT INTEGER(I-N) 19 : -INC CCOPTIO 20 : -INC SMELEME 21 : -INC SMRIGID 22 : -INC SMMATRI 23 : -INC CCREEL 24 : SEGMENT,INUINV(NNGLOB) 25 : SEGMENT,ITOPO(IENNO) 26 : SEGMENT,IITOP(NNOE+1) 27 : SEGMENT,IMINI(INC) 28 : SEGMENT,IPOS(NNOE1) 29 : SEGMENT,INCTRR(NIRI) 30 : SEGMENT,INCTRA(NLIGRE) 31 : SEGMENT,IPV(NNOE) 32 : C 33 : C **** CES TABLEAUX SERVENT AU REPERAGE DE LA MATRICE POUR L'ASSEMBLAG 34 : C **** IL SERONT TOUS SUPPRIMES EN FIN D'ASSEMBLAGE. 35 : C 36 : ** 37 : SEGMENT,IVAL(NNN) 38 : SEGMENT,ITRA(NNN,2) 39 : SEGMENT TRATRA 40 : REAL*8 XTRA(INCRED,INCDIF) 41 : INTEGER LTRA(INC,INCDIF) 42 : INTEGER NTRA(INCRED,INCDIF) 43 : INTEGER MTRA(INCDIF) 44 : ENDSEGMENT 45 : C **** IVAL(I)=J : LA I EME LIGNE D'UNE PETITE MATRICE S'ASSEMBLE 46 : C DANS LA J EME DE LA GRANDE. 47 : C **** ITRAV(I,1)=J : LA IEME INCONNUE DU NOEUD EN COURS D'ASSEMBLAGE 48 : C ET QUI SE TROUVE DANS LA PETITE MATRICE SE TROUVE 49 : C EN J EME POSITION DE LA PETITE MATRICE. 50 : C **** ITRAV(I,2) : LA IEME INCONNUE DU NOEUD EN COURS D'ASSEMBLAGE 51 : C PRESENT DANS LA PETITE MATRICE EST EN JEME 52 : C POSITION DANS LA GRANDE 53 : SEGMENT,RA(N1,N1)*D 54 : SEGMENT JNOMUL 55 : LOGICAL INOMUL(NNR) 56 : ENDSEGMENT 57 : REAL*8 DMAX,COER,DDDD,DMAXY,DMAXGE 58 : LOGICAL NOMUL 59 : * 60 : * PV ON ACTIVE UNE FOIS POUR TOUTES LES MELEME DESCR... DE LA RIGIDITE 61 : * ON EN PROFITE POUR CREER INOMUL 62 : C 63 : C **** RECHERCHE DE LA DIMENSION MAX DE IVAL,ET SEGINI DE IVAL ET ITRA 64 : C 65 : INCTRR=INCTR1 66 : SEGACT,INCTRR 67 : MRIGID=ITRAV1 68 : SEGACT,MRIGID 69 : NNR=IRIGEL(/2) 70 : NNN=0 71 : SEGINI JNOMUL 72 : DO 1 IRI=1,NNR 73 : INCTRA=INCTRR(IRI) 74 : SEGACT INCTRA 75 : DESCR=IRIGEL(3,IRI) 76 : SEGACT, DESCR 77 : IPT1=IRIGEL(1,IRI) 78 : SEGACT IPT1 79 : * XMATRI=IRIGEL(4,IRI) 80 : * SEGACT XMATRI 81 : NA=LISINC(/2) 82 : NNN=MAX(NA,NNN) 83 : INOMUL(IRI)=.TRUE. 84 : IF(IPT1.ITYPEL.EQ.22) INOMUL(IRI)=.FALSE. 85 : 1 CONTINUE 86 : SEGINI,IVAL 87 : SEGINI,ITRA 88 : C 89 : C **** ACTIVATION DES SEGMENTS DE TRAVAILS ET DE MMATRI 90 : C 91 : IMINI=IMINI1 92 : SEGACT,IMINI 93 : N1=IMINI(/1) 94 : ITOPO=ITOPO1 95 : SEGACT,ITOPO 96 : IITOP=IITOP1 97 : SEGACT,IITOP 98 : INUINV=INUIN1 99 : SEGACT,INUINV 100 : IPOS=IPO1 101 : SEGACT,IPOS 102 : NNOE=IPOS(/1)-1 103 : INC=IPOS(NNOE+1) 104 : MMATRI=MMMTRI 105 : SEGACT,MMATRI*MOD 106 : SEGINI,MDIAG 107 : SEGINI,MILIGN 108 : IILIGN=MILIGN 109 : IDIAG=MDIAG 110 : MINCPO=IINCPO 111 : SEGACT,MINCPO 112 : INCDIF=INCPO(/1) 113 : MIMIK=IIMIK 114 : SEGACT,MIMIK 115 : SEGINI IPV 116 : INCRED=0 117 : DO 80 INO =1,NNOE 118 : ICOMPT=0 119 : MAXELE = (IITOP(INO+1)-IITOP(INO))/2 120 : DO 81 IELE=1,MAXELE 121 : IIU=IITOP(INO) + IELE + IELE -2 122 : IEL=ITOPO(IIU) 123 : IRI=ITOPO(IIU+1) 124 : MELEME=IRIGEL(1,IRI) 125 : DO 83 I=1,NUM(/1) 126 : IP=INUINV(NUM(I,IEL)) 127 : IF (IP.GT.INO) GOTO 83 128 : IF (IPV(IP).EQ.INO) GOTO 83 129 : IPV(IP)=INO 130 : ICOMPT=ICOMPT+1 131 : 83 CONTINUE 132 : 81 CONTINUE 133 : INCRED=MAX(INCRED,ICOMPT) 134 : 80 CONTINUE 135 : SEGSUP IPV 136 : INCRED=INCRED*INCDIF 137 : SEGINI TRATRA 138 : C 139 : C **** BOUCLE *100* SUR LES NUMEROS DE NOEUDS QUE L'ON ASSEMBLE 140 : C 141 : LLVNUL=0 142 : IJMAX=0 143 : NJTOT=0 144 : NNOE=IPOS(/1)-1 145 : DNORME=0. 146 : SEGINI MDNOR 147 : IDNORM=MDNOR 148 : * 149 : * en cas de normalisation des variables. 150 : * 151 : IF(NORINC.NE.0) THEN 152 : inwuit=0 152 : CALL ASSE10(ITRAV1,1,MDNOR,MIMIK,MINCPO,INUIN1,inwuit) 153 : ELSE 154 : DO 53 IU=1,INC 155 : DNOR(IU)=1.d0 156 : 53 CONTINUE 157 : ENDIF 158 : * On balaye les noeuds dans l'ordre des elements 159 : DO 100 JR=1,IRIGEL(/2) 160 : IPT6=IRIGEL(1,JR) 161 : DO 100 JL=1,IPT6.num(/2) 162 : DO 100 JP=1,IPT6.num(/1) 163 : INO=INUINV(IPT6.NUM(JP,JL)) 164 : ** DO 100 INO=1,NNOE 165 : IF (ILIGN(INO).NE.0) GOTO 100 166 : DO 101 IIT=1,INCDIF 167 : MTRA(IIT)=0 168 : 101 CONTINUE 169 : IPRE=IPOS(INO)+1 170 : IDER=IPOS(INO+1) 171 : LLVVA=0 172 : C 173 : C **** BOUCLE *99* SUR LES ELEMENTS TOUCHANT LE NOEUD INO 174 : C POUR LES ELEMNTS MULTIPLICATEUR ON NE FAIT PAS 175 : C L'ASSEMBLAGE 176 : C 177 : MAXELE= (IITOP(INO+1) -IITOP(INO))/2 178 : DO 99 IELE=1,MAXELE 179 : IIU=IITOP(INO) + IELE + IELE - 2 180 : IEL=ITOPO(IIU) 181 : IRI=ITOPO(IIU+1) 182 : MELEME=IRIGEL(1,IRI) 183 : DESCR=IRIGEL(3,IRI) 184 : INCTRA=INCTRR(IRI) 185 : XMATRI=IRIGEL(4,IRI) 186 : SEGACT XMATRI 187 : COER=COERIG(IRI) 188 : C 189 : C **** NOMUL =.FALSE. IL EXISTE UN MULTUIPLICATEUR 190 : C **** INITIALISATION DE IVAL. IVAL(I)=J VEUT DIRE QUE 191 : C **** LA I EME LIGNE DE LA PETITE MATRICE S'ASSEMBLE DANS 192 : C **** LA J EME DE LA GRANDE MATRICE. 193 : C 194 : NIN=LISINC(/2) 195 : NOMUL=INOMUL(IRI) 196 : NA=0 197 : DO 98 ICO=1,NIN 198 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 199 : IJB=INCTRA(ICO) 200 : IVAL(ICO)=INCPO(IJB,IJA) 201 : IF(IJA.NE.INO) GO TO 98 202 : NA=NA+1 203 : ITRA(NA,1)=ICO 204 : ITRA(NA,2)=IVAL(ICO) 205 : 98 CONTINUE 206 : * XMATRI=IMATTT(IEL) 207 : * SEGACT,XMATRI 208 : C 209 : C **** BOUCLE *85* SUR LES INCONNUES DE LA PETITE MATRICE 210 : C 211 : DO 95 IK=1,NIN 212 : IO=IVAL(IK) 213 : IF(IO.GT.IDER) GO TO 95 214 : IPOO=IK*NIN - NIN 215 : DO 90 INCC=1,NA 216 : INCO=ITRA(INCC,2) 217 : IF(IO.GT.INCO) GO TO 90 218 : ILOC=INCO-IPRE+1 219 : JJ=ITRA(INCC,1) 220 : * IF(IK.LE.JJ) THEN 221 : * IPO=(JJ*(JJ-1))/2+IK 222 : * ELSE 223 : * IPO=(IK*(IK-1))/2+JJ 224 : * ENDIF 225 : IPO=IPOO+JJ 226 : ILTT= LTRA(IO,ILOC) 227 : IF(ILTT.EQ.0) THEN 228 : LLVVA=LLVVA+1 229 : IMMTT=MTRA(ILOC)+1 230 : MTRA(ILOC)=IMMTT 231 : XTRA(IMMTT,ILOC)=0.D0 232 : NTRA(IMMTT,ILOC)=IO 233 : LTRA(IO,ILOC)=IMMTT 234 : ILTT=IMMTT 235 : ENDIF 236 : IF(NOMUL) THEN 237 : XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IK,JJ,IEL)*COER 238 : ENDIF 239 : 90 CONTINUE 240 : 95 CONTINUE 241 : SEGDES,XMATRI 242 : 99 CONTINUE 243 : C 244 : C *** COMPACTAGE DES LIGNES, EN MEME TEMPS CALCUL DE IJMAXQUI SERA 245 : C *** LA DIMENSION MAX D'UN SEGMENT LIGN. 246 : C *** LE SEGMENT ASSOCIE A UNE LIGNE (SEGMENT LLIGN)EST DE LA FORME : 247 : C *** IMMMM(NA) PERMET DE SAVOIR SI UN MOUVENENT D'ENSEMBLE SUR LA 248 : C *** LIGNE EXISTE. IPPO(NA+1) DONNE LA POSITION DANS XXVA LA 1ERE 249 : C *** VALEUR DE LA LIGNE .XXVA VALEUR DE LA MATRICE. 250 : C *** LINC(I)DONNE LE NUMERO DE LA COLONNE DU IEME ELEM DE XXVA 251 : C 252 : NA = IDER-IPRE+1 253 : LLVNUL=LLVNUL+LLVVA 254 : SEGINI,LLIGN 255 : ILIGN(INO)=LLIGN 256 : NBA=0 257 : DO 120 JPA=1,NA 258 : IIIN=IPRE+JPA -1 259 : IMMMM(JPA)=IIIN 260 : IPPO(JPA)=NBA 261 : DO 121 IPAK = 1,MTRA(JPA) 262 : IUNPAK=NTRA(IPAK,JPA) 263 : LTRA(IUNPAK,JPA)=0 264 : NBA=NBA+1 265 : LINC(NBA)=IUNPAK 266 : XXVA(NBA)=XTRA(IPAK,JPA) 267 : IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA) 268 : 121 CONTINUE 269 : 120 CONTINUE 270 : IPPO(NA+1)= NBA 271 : NJMAX= 0 272 : DO 126 JPA=IPRE,IDER 273 : IPNO(JPA)=INO 274 : IPDE=IPPO(JPA-IPRE+1)+1 275 : IPDF=IPPO(JPA-IPRE+2) 276 : LPA=JPA 277 : DO 155 JHT=IPDE,IPDF 278 : LPA=MIN(LPA,LINC(JHT)) 279 : 155 CONTINUE 280 : LDEB(JPA-IPRE+1)=LPA 281 : NNA= JPA- LPA +1 282 : NJMAX=NJMAX+NNA 283 : 126 CONTINUE 284 : NJTOT=NJTOT+NJMAX 285 : IF(IJMAX.LT.NJMAX) IJMAX=NJMAX 286 : SEGDES,LLIGN 287 : 100 CONTINUE 288 : SEGSUP TRATRA 289 : C 290 : C **** ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS 291 : C **** POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE 292 : C **** A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE. 293 : C 294 : * d'abord etablir une norme generale pour le cas ou on n'arrive pas 295 : * a calculer la norme particuliere 296 : DMAXGE=1D-31 297 : DO 378 I=1,INC 298 : DMAXGE=MAX(DMAXGE,ABS(DIAG(I))) 299 : 378 CONTINUE 300 : IF(DMAXGE.LE.1D-30) DMAXGE = 1.D0 301 : if (iimpi.eq.1524) 302 : > write (6,*) ' nb inconnues facteur multiplicatif general ', 303 : > INC,DMAXGE 304 : IENMU=0 305 : 375 IENMU1 = IENMU 306 : IENMU=0 307 : DO 376 I=1,NNR 308 : IF(.NOT.INOMUL(I)) IENMU=IENMU+1 309 : 376 CONTINUE 310 : IF( IENMU.EQ.0) GO TO 3750 311 : MIMIK=IIMIK 312 : SEGACT,MIMIK 313 : DO 11 I=1,NNR 314 : IF(INOMUL(I)) GO TO 11 315 : DESCR=IRIGEL(3,I) 316 : N3=LISINC(/2) 317 : COER=COERIG(I) 318 : MELEME=IRIGEL(1,I) 319 : INCTRA=INCTRR(I) 320 : XMATRI=IRIGEL(4,I) 321 : SEGACT XMATRI 322 : N2=NUM(/2) 323 : IF (RE(/3).EQ.0) THEN 324 : INOMUL(I)=.TRUE. 325 : SEGDES XMATRI 326 : GOTO 11 327 : ENDIF 328 : * XMATRI=IMATTT(1) 329 : * SEGACT,XMATRI 330 : N1=RE(/1) 331 : SEGINI,RA 332 : DO 14 IEL=1,N2 333 : DO 15 ICO=1,N3 334 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 335 : IJB=INCTRA(ICO) 336 : IVAL(ICO)=INCPO(IJB,IJA) 337 : 15 CONTINUE 338 : DDDD= 1D-31 339 : DMAX=DDDD 340 : if (n3.lt.3) then 341 : do ico=1,n3 342 : DMAX=MAX(DMAX,ABS(DIAG(IVAL(ICO)))) 343 : enddo 344 : endif 345 : DO 19 ICO=3,N3 346 : DMAX=MAX(DMAX,ABS(DIAG(IVAL(ICO)))) 347 : 19 CONTINUE 348 : * write (6,*) ' assem2 dmax dmaxge ',dmax,dmaxge 349 : C AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION 350 : IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7391)DMAX,DDDD,IENMU,IENMU1 351 : 1,I,IEL 352 : 7391 FORMAT(' DMAX DDDD IENMU IENMU1 I IEL',2E12.5,4I3) 353 : IF(DMAX.LE.1D-30.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1) GO TO 377 354 : IF(DMAX.LE.1D-30) DMAX = DMAXGE 355 : DMAX=DMAX*1.5 356 : * on penalise la matrice en cas de resolution iterative 357 : **pv if (nucrou.eq.1) DMAX=DMAX*1D5 358 : ** write (6,*) ' assem2 i iel dmax ',i,iel,dmax 359 : * XMATRI=IMATTT(IEL) 360 : * SEGACT,XMATRI 361 : DMAXY=SQRT(XPETIT)*1D5 362 : if (norinc.eq.0) dmaxy=1.D0 363 : DO 821 ICO=3,N1 364 : * DO 82110 IKO=1,N1 365 : DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,IEL))) 366 : *82110 CONTINUE 367 : 821 CONTINUE 368 : * if (dmaxy.lt.1d-50) write (6,*) (re(ico,1),ico=1,n1) 369 : * if (dmaxy.lt.1d-50)write (6,*)' assem2 dmax dmaxy ',dmax,dmaxy,n1 370 : DMAX = DMAX / DMAXY 371 : IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7398) DMAX 372 : 7398 FORMAT(' facteur multiplicatif de norme ',e12.5) 373 : DO 21 ICO=1,N1 374 : DO 2110 IKO=1,N1 375 : RA(ICO,IKO)=RE(ICO,IKO,IEL)*DMAX*COER 376 : 2110 CONTINUE 377 : 21 CONTINUE 378 : RA(1,1)=RA(1,1)*DMAXY 379 : RA(2,1)=RA(2,1)*DMAXY 380 : RA(1,2)=RA(1,2)*DMAXY 381 : RA(2,2)=RA(2,2)*DMAXY 382 : DO 22 ICO=1,2 383 : DNOR(IVAL(ICO))=DMAX 384 : 22 CONTINUE 385 : if (abs(dmax).gt. 1d50) write (6,*) ' assem2 dmax dmaxy ',dmax, 386 : > dmaxy,ival(ico) 387 : DO 24 ICO=1,N3 388 : INO=INUINV(NUM(NOELEP(ICO),IEL)) 389 : IO=IVAL(ICO) 390 : if (ico.eq.1) io1=io 391 : if (ico.eq.2) io2=io 392 : LLIGN=ILIGN(INO) 393 : SEGACT,LLIGN*MOD 394 : * DIAG(IO)=DIAG(IO)+RA((ICO*(ICO+1))/2) 395 : DIAG(IO)=DIAG(IO)+RA(ICO,ICO) 396 : DO 132 JLIJ=1,IMMMM(/1) 397 : JLIJ1=JLIJ 398 : IF( IMMMM(JLIJ).EQ.IO) GO TO 133 399 : 132 CONTINUE 400 : IF(IIMPI.EQ.1524) WRITE(IOIMP,7354) 401 : 7354 FORMAT( ' PREMIERE ERREUR 5') 401 : CALL ERREUR(5) 402 : RETURN 403 : 133 CONTINUE 404 : * IREE=ICO*N3-ICO 405 : DO 26 IRO=1,N3 406 : IA=IVAL(IRO) 407 : IF(IA.GT.IO) GO TO 26 408 : * IF(IRO.LE.ICO) IRE=(ICO*(ICO-1))/2+IRO 409 : * IF(IRO.GT.ICO) IRE=(IRO*(IRO-1))/2+ICO 410 : JLT=IPPO(JLIJ1+1) 411 : JLD=IPPO(JLIJ1)+1 412 : DO 134 JL=JLD,JLT 413 : JL1=JL 414 : IF(LINC(JL).EQ.IA) GO TO 135 415 : 134 CONTINUE 416 : IF(IIMPI.NE.1524) WRITE(IOIMP,7355) 417 : 7355 FORMAT( ' DEUXIEME ERREUR 5') 417 : CALL ERREUR(5) 418 : RETURN 419 : 135 CONTINUE 420 : XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO) 421 : 26 CONTINUE 422 : SEGDES,LLIGN 423 : 24 CONTINUE 424 : * on stocke dans ittr les couples de LX 425 : ittr(io1)=io2 426 : ittr(io2)=io1 427 : * SEGDES,XMATRI 428 : 14 CONTINUE 429 : INOMUL(I)=.TRUE. 430 : 377 CONTINUE 431 : SEGSUP,RA 432 : 11 CONTINUE 433 : GO TO 375 434 : 3750 CONTINUE 435 : DO 18 IK=1,NNR 436 : INCTRA=INCTRR(IK) 437 : SEGSUP,INCTRA 438 : 18 CONTINUE 439 : * PV ON DESACTIVE TOUT 440 : NNR=IRIGEL(/2) 441 : DO 2 IRI=1,NNR 442 : DESCR=IRIGEL(3,IRI) 443 : SEGDES DESCR 444 : IPT1=IRIGEL(1,IRI) 445 : SEGDES IPT1 446 : XMATRI=IRIGEL(4,IRI) 447 : SEGDES XMATRI 448 : 2 CONTINUE 449 : 450 : INTERR(1)=NJTOT 451 : IF(IIMPI.EQ.1457) WRITE(IOIMP,4821) LLVNUL,NJTOT 452 : 4821 FORMAT(' NB DE VALEURS NON NULLES DANS LA MATRICE ',I9,/ 453 : # ' NB DE VALEURS DANS LA MATRICE ',I9) 454 : IF(NORINC.NE.0) THEN 454 : CALL ASSE10(MRIGID,2,MDNOR,MIMIK,MINCPO,INUIN1,inwuit) 455 : SEGDES,MRIGID 456 : ELSE 457 : SEGDES,MRIGID 458 : ENDIF 459 : SEGSUP,INCTRR 460 : SEGDES,MDIAG 461 : SEGDES,MIMIK 462 : SEGDES,MDNOR 463 : SEGSUP,IMINI 464 : SEGSUP,ITOPO 465 : SEGSUP,IITOP 466 : SEGDES,MILIGN 467 : SEGSUP,INUINV 468 : SEGDES,MMATRI 469 : MMMTRI=MMATRI 470 : SEGDES,MINCPO 471 : SEGDES,IPOS 472 : SEGSUP,IVAL,ITRA,JNOMUL 473 : RETURN 474 : END 475 : 476 : 477 : 478 : 479 : 480 : 481 : 482 : 483 : 484 : 485 : 486 : 487 : 488 : 489 : 490 : 491 : 492 : 493 : 494 : 495 : 496 :
© Cast3M 2003 - All rights reserved.
Disclaimer