1 : C ASSEM5 SOURCE PV 11/05/03 21:15:16 6959 2 : c assem5 source chat 94/08/12 21:15:05 1247 3 : SUBROUTINE ASSEM5(ITRAV1,ITOPO1,INUIN1,MMMTRI,INCTR1 4 : &,IITOP1,NBNNMA,SNTT) 5 : c 6 : cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 7 : c 8 : c **** subroutine pour faire l'assemblage de matrices symetriques 9 : c en vue d'une decomposition de crout modifiee. 10 : c 11 : c en entree: 12 : c **** itrav1 : pointeur objet mrigidite 13 : c **** itopo1 : pointeur segment de travail itopo ( voir assem1) 14 : c **** iitop1 : pointeur segment de travail iitop ( voir assem1) 15 : c **** inuin1 : pointeur segment de travail inuinv( voir assem1) 16 : c **** mmmtri : pointeur objet matrice triangularisee (non modifie) 17 : c (voir smmatri) 18 : c **** nbnnma : nombre d'inconnues maitresses 19 : c **** sntt : segment des noeuds non totalement maitres 20 : c 21 : cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 22 : c 23 : IMPLICIT REAL*8 (A-H,O-Z) 24 : IMPLICIT INTEGER(I-N) 25 : -INC CCREEL 26 : -INC CCOPTIO 27 : -INC SMELEME 28 : -INC SMRIGID 29 : -INC SMMATRI 30 : c 31 : SEGMENT,INUINV(NNGLOB) 32 : SEGMENT,ITOPO(IENNO) 33 : SEGMENT,IITOP(NNOE+1) 34 : SEGMENT,INCTRR(NIRI) 35 : SEGMENT,INCTRA(NLIGRE) 36 : SEGMENT,IPV(NNOE) 37 : SEGMENT SNTT 38 : INTEGER NTTMAI(NN) 39 : ENDSEGMENT 40 : c 41 : c **** ces tableaux servent au reperage de la matrice pour l'assemblag 42 : c **** il seront tous supprimes en fin d'assemblage. 43 : c 44 : SEGMENT,IVAL(NNN) 45 : SEGMENT,ITRA(NNN,2) 46 : SEGMENT TRATRA 47 : REAL*8 XTRA(INCRED,INCDIF) 48 : INTEGER LTRA(INC,INCDIF) 49 : INTEGER NTRA(INCRED,INCDIF) 50 : INTEGER MTRA(INCDIF) 51 : ENDSEGMENT 52 : c 53 : cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 54 : c 55 : c **** ival(i)=j : la i eme ligne d'une petite matrice s'assemble 56 : c dans la j eme de la grande. 57 : c **** itrav(i,1)=j : la ieme inconnue du noeud en cours d'assemblage 58 : c et qui se trouve dans la petite matrice se trouve 59 : c en j eme position de la petite matrice. 60 : c **** itrav(i,2) : la ieme inconnue du noeud en cours d'assemblage 61 : c present dans la petite matrice est en jeme 62 : c position dans la grande 63 : c 64 : cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 65 : c 66 : SEGMENT,RA(N1,N1)*D 67 : SEGMENT JNOMUL 68 : LOGICAL INOMUL(NNR) 69 : ENDSEGMENT 70 : REAL*8 DMAX,COER,DDDD,DMAXY 71 : LOGICAL NOMUL 72 : * 73 : * en cas de normalisation des variables. 74 : * 75 : INWUIT=0 75 : IF(NORINC.NE.0) CALL ASSEM0(ITRAV1,1,INWUIT) 76 : * 77 : * pv on active une fois pour toutes les meleme descr... de la rigidite 78 : * on en profite pour creer inomul 79 : c 80 : c **** recherche de la dimension max de ival,et segini de ival et itra 81 : c 82 : c 83 : INCTRR=INCTR1 84 : SEGACT,INCTRR 85 : c 86 : MRIGID=ITRAV1 87 : 88 : SEGACT,MRIGID 89 : NNR=IRIGEL(/2) 90 : NNN=0 91 : SEGINI JNOMUL 92 : DO 1 IRI=1,NNR 93 : c 94 : INCTRA=INCTRR(IRI) 95 : SEGACT INCTRA 96 : c 97 : DESCR=IRIGEL(3,IRI) 98 : SEGACT, DESCR 99 : c 100 : IPT1=IRIGEL(1,IRI) 101 : SEGACT IPT1 102 : c 103 : XMATRI=IRIGEL(4,IRI) 104 : SEGACT XMATRI 105 : c 106 : c ... na = nb de variables primales ... 107 : NA=LISINC(/2) 108 : c ... nnn = max des nb de variables primales ... 109 : NNN=MAX(NA,NNN) 110 : c 111 : c ... inomul (non multiplicateurs ?) dit s'il ne s'agit pas des cl ... 112 : INOMUL(IRI)=.TRUE. 113 : IF(IPT1.ITYPEL.EQ.22) INOMUL(IRI)=.FALSE. 114 : c 115 : 1 CONTINUE 116 : c 117 : c ... ival a la taille nnn ... 118 : SEGINI,IVAL 119 : c ... itra a la taille (nnn,2) ... 120 : SEGINI,ITRA 121 : c 122 : c **** activation des segments de travails et de mmatri 123 : c 124 : c ... recherche du nombre d'inconnues impliquees 125 : c 126 : MMATRI=MMMTRI 127 : SEGACT,MMATRI*MOD 128 : c 129 : MINCPO=IINCPO 130 : SEGACT,MINCPO 131 : c 132 : c 133 : MAXIPO=0 134 : DO 150 I=1,INCPO(/2) 135 : DO 150 J=1,INCPO(/1) 136 : MAXIPO=MAX(INCPO(J,I),MAXIPO) 137 : 150 CONTINUE 138 : c 139 : N1=MAXIPO 140 : c 141 : ITOPO=ITOPO1 142 : SEGACT,ITOPO 143 : c 144 : IITOP=IITOP1 145 : SEGACT,IITOP 146 : c 147 : INUINV=INUIN1 148 : SEGACT,INUINV 149 : c 150 : c ... nnoe = nombre de noeuds impliqués ... 151 : c 152 : SEGACT SNTT 153 : NNOE=INCPO(/2)+NTTMAI(/1) 154 : c 155 : c ... inc = dimension de la matrice ... 156 : INC=MAXIPO 157 : c 158 : SEGINI,MDIAG 159 : IDIAG=MDIAG 160 : c 161 : SEGINI,MILIGN 162 : IILIGN=MILIGN 163 : c 164 : c ... incdif = nb de variables primales différentes ... 165 : INCDIF=INCPO(/1) 166 : c 167 : MIMIK=IIMIK 168 : SEGACT,MIMIK 169 : c 170 : NNOE=INCPO(/2) 171 : c 172 : c ... ipv est de taille nnoe ... 173 : SEGINI IPV 174 : c 175 : INCRED=0 176 : DO 80 INO =1,NNOE 177 : ICOMPT=0 178 : c ... maxele = nb d'éléments qui contiennent le noeud no ino ... 179 : MAXELE = (IITOP(INO+1)-IITOP(INO))/2 180 : c ... boucle sur ces éléments ... 181 : DO 81 IELE=1,MAXELE 182 : c 183 : c ... iiu = pointeur dans iitop ... 184 : IIU=IITOP(INO) + IELE + IELE -2 185 : c 186 : c ... iel = numero de l'élément dans son maillage ... 187 : IEL=ITOPO(IIU) 188 : c 189 : c ... iri = numero du maillage (dans irigel) de cet élément ... 190 : IRI=ITOPO(IIU+1) 191 : c 192 : c ... on va chercher ce maillage dans irigel ... 193 : MELEME=IRIGEL(1,IRI) 194 : c 195 : c ... boucle sur les noeuds de cet élément ... 196 : DO 83 I=1,NUM(/1) 197 : c 198 : c ... ip = numéro local du i-ème noeud de l'élément ... 199 : IP=INUINV(NUM(I,IEL)) 200 : c 201 : c ... la première condition = les "connections" sont symétriques ... 202 : **pv IF (IP.GT.INO) GOTO 83 203 : c ... le noeud ip a déjà vu ce noeud ... 204 : IF (IPV(IP).EQ.INO) GOTO 83 205 : c 206 : c ... ces deux opérations se font si (ip <= ino et ipv(ip) != ino) ... 207 : c ... ipv = n° max des noeuds "connectés" au noeud ip ... 208 : c 209 : IPV(IP)=INO 210 : c ... icompt = nombre de noeuds connectés (de numéro < ou =) au noeud ino ... 211 : ICOMPT=ICOMPT+1 212 : 83 CONTINUE 213 : 81 CONTINUE 214 : c ... incred = maximum des icompt sur tous les noeuds impliqués ... 215 : c 216 : INCRED=MAX(INCRED,ICOMPT) 217 : 80 CONTINUE 218 : c 219 : SEGSUP IPV 220 : c 221 : INCRED=INCRED*INCDIF 222 : c ... la taille de tratra dépend de inc, incred et incdif ... 223 : SEGINI TRATRA 224 : c 225 : c ... coefficients de normalisation ... 226 : c ----------------------------- 227 : c 228 : LLVNUL=0 229 : IJMAX=0 230 : NJTOT=0 231 : c 232 : c ... la taille de mdnor est inc (nombre d'équations) ... 233 : SEGINI MDNOR 234 : IDNORM=MDNOR 235 : DO 20 I=1,INC 236 : DNOR(I)=1.D0 237 : 20 CONTINUE 238 : c 239 : c ... si une normalisation est imposée, les valeurs dans mdnor seront modifiées 240 : c par assem3 ... 240 : IF(NORINC.NE.0) CALL ASSEM3(MDNOR,MIMIK,MINCPO) 241 : c 242 : c **** boucle *100* sur les numeros de noeuds que l'on assemble 243 : c ------------------------------------------------------------- 244 : c 245 : c ... ipre = n° de la première ligne concernée par le noeud ino ... 246 : IPRE=1 247 : c 248 : DO 999 IJK=1,2 249 : c 250 : DO 100 INO=1,NNOE 251 : DO 101 IIT=1,INCDIF 252 : MTRA(IIT)=0 253 : 101 CONTINUE 254 : c ... ider = n° de la dernière ligne concernée par le noeud ino ... 255 : IDER=0 256 : c 257 : IF(IJK.EQ.2) IPRE=MAXIPO+1 258 : DO 998 II=1,INCDIF 259 : IF (IJK.EQ.1.AND.INCPO(II,INO).LE.NBNNMA) 260 : # IDER=MAX(INCPO(II,INO),IDER) 261 : IF (IJK.EQ.2.AND.INCPO(II,INO).GT.NBNNMA) THEN 262 : IDER=MAX(INCPO(II,INO),IDER) 263 : IPRE=MIN(INCPO(II,INO),IPRE) 264 : ENDIF 265 : 998 CONTINUE 266 : IF (IDER.EQ.0.OR.IPRE.EQ.(MAXIPO+1)) GOTO 100 267 : LLVVA=0 268 : c 269 : c **** boucle *99* sur les elements touchant le noeud ino 270 : c pour les elements multiplicateur on ne fait pas 271 : c l'assemblage 272 : c 273 : c ... maxele = nb d'éléments qui contiennent le noeud no ino ... 274 : MAXELE= (IITOP(INO+1) -IITOP(INO))/2 275 : c 276 : 277 : DO 99 IELE=1,MAXELE 278 : IIU=IITOP(INO) + IELE + IELE - 2 279 : c ... iel = numero de l'élément dans son maillage ... 280 : IEL=ITOPO(IIU) 281 : c ... iri = numero du maillage (dans irigel) de cet élément ... 282 : IRI=ITOPO(IIU+1) 283 : c ... meleme = pointeur vers ce maillage ... 284 : MELEME=IRIGEL(1,IRI) 285 : c ... descr = pointeur vers le segment descr concerné ... 286 : DESCR=IRIGEL(3,IRI) 287 : c ... inctra contient les indices (dans imik et ihar) des descriptions 288 : c des ddl correspondants ... 289 : INCTRA=INCTRR(IRI) 290 : c ... xmatri = tableau de pointeurs vers les matrices élémentaires ... 291 : XMATRI=IRIGEL(4,IRI) 292 : 293 : c ... coer = coefficient multiplicateur ... 294 : COER=COERIG(IRI) 295 : c 296 : c **** nomul =.false. il existe un multiplicateur 297 : c **** initialisation de ival. ival(i)=j veut dire que 298 : c **** la i eme ligne de la petite matrice s'assemble dans 299 : c **** la j eme de la grande matrice. 300 : c 301 : c ... nin = nombre de variables primales (dans descr) ... 302 : NIN=LISINC(/2) 303 : c ... non multiplicateur ? ... 304 : NOMUL=INOMUL(IRI) 305 : NA=0 306 : c ... boucle sur les variables primales de l'élément ... 307 : c 308 : DO 98 ICO=1,NIN 309 : c ... ija = numéro local du noeud-support du ddl ... 310 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 311 : c ... ijb = indice du ddl ... 312 : IJB=INCTRA(ICO) 313 : c ... on met dans ival(ico) le n° de l'équation correspondant au noeud ija 314 : c et ddl no ijb ... 315 : IVAL(ICO)=INCPO(IJB,IJA) 316 : c ... si on a trouvé le bon noeud ... 317 : IF(IJA.EQ.INO) THEN 318 : c ... on incrémente na (le nombre de ddl connectés au noeud ino ?) ... 319 : NA=NA+1 320 : ITRA(NA,1)=ICO 321 : ITRA(NA,2)=IVAL(ICO) 322 : ENDIF 323 : 98 CONTINUE 324 : * XMATRI=IMATTT(IEL) 325 : * SEGACT,XMATRI 326 : c 327 : c **** boucle *85* sur les inconnues de la petite matrice 328 : c 329 : DO 95 IK=1,NIN 330 : c ... io = n° de l'équation correspondant au ddl no ik de l'élément (n° de la colonne) ... 331 : IO=IVAL(IK) 332 : c ... s'il dépasse le dernier intéressant on ne fait rien ... 333 : IF(IO.GT.IDER) GO TO 95 334 : c ... boucle sur les ddl's du noeud ino ... 335 : DO 90 INCC=1,NA 336 : c ... inco = le n° d'équation (ou de ligne) du ddl no incc du noeud ino ... 337 : INCO=ITRA(INCC,2) 338 : c 339 : IF (IJK.EQ.1) THEN 340 : IF (INCO.GT.IDER) GOTO 90 341 : ENDIF 342 : c 343 : c ... si io > inco on ne fait plus rien dans cette boucle ... 344 : c ... car c'est le cas symétrique !!! ... 345 : IF (IO.GT.INCO) GO TO 90 346 : c 347 : c ... iloc = n° de ligne relatif / ipre du noeud ino ... 348 : c 349 : ILOC=INCO-IPRE+1 350 : IF (IJK.EQ.2) THEN 351 : IF (ILOC.LE.0) GOTO 90 352 : ENDIF 353 : c 354 : c ... jj = n° (dans l'élément) de la colonne (et de la ligne) correspondant au 355 : c ddl n° incc du noeud ino ... 356 : JJ=ITRA(INCC,1) 357 : c ... ci-dessous on stockera les quelques lignes de la rigidité concernant le 358 : c noeud ino dans le segment tratra. les coefficients vont dans xtra, le deuxième 359 : c indice (iloc) est le numéro relatif de la ligne, le premier (imtt) n'a rien à voir avec le 360 : c numéro de la colonne. celui-ci est stocké dans ntra(imtt,iloc). l'information inverse 361 : c (imtt en fonction de io) se trouve dans ltra(io,iloc). on met dans mtra le nombre 362 : c de termes différents assemblés par ligne ... 363 : ILTT= LTRA(IO,ILOC) 364 : IF(ILTT.EQ.0) THEN 365 : c ... llvva est initialisé à 0 au début de la boucle 100 ... 366 : LLVVA=LLVVA+1 367 : c ... les mtra sont mis à 0 au début de la boucle 100 ... 368 : IMMTT=MTRA(ILOC)+1 369 : MTRA(ILOC)=IMMTT 370 : XTRA(IMMTT,ILOC)=0.D0 371 : NTRA(IMMTT,ILOC)=IO 372 : LTRA(IO,ILOC)=IMMTT 373 : ILTT=IMMTT 374 : ENDIF 375 : c 376 : c ... si ce n'est pas une condition aux limites, on assemble ... 377 : IF(NOMUL) THEN 378 : cfaux !!! xtra(iltt,iloc)=xtra(iltt,iloc)+re(ik,jj)*coer 379 : c ... cette ligne a été remplacée par ce qui suit car pour les termes 380 : c de couplage entre les ddl du même noeud, les termes au-dessus de 381 : c la diagonale étaient pris. Ça marchait car les matrices étaient 382 : c symétriques, comme maintenant on travaillera aussi sur des asymétriques ... 383 : IF(JJ.LT.IK) THEN 384 : IIILIG=IK 385 : IIICOL=JJ 386 : ELSE 387 : IIILIG=JJ 388 : IIICOL=IK 389 : ENDIF 390 : c ... attention ! dans xtra 1er indice est lié à la colonne, deuxième est le numéro 391 : c de la ligne, tandis que dans re c'est l'inverse ... 392 : XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IIILIG,IIICOL,IEL) 393 : & *COER 394 : ENDIF 395 : 90 CONTINUE 396 : 95 CONTINUE 397 : * SEGDES,XMATRI 398 : 99 CONTINUE 399 : c 400 : c *** compactage des lignes, en meme temps calcul de ijmax qui sera 401 : c *** la dimension max d'un segment lign. 402 : c *** le segment associe a une ligne (segment llign) est de la forme : 403 : c *** immmm(na) permet de savoir si un mouvement d'ensemble sur la 404 : c *** ligne existe. ippo(na+1) donne la position dans xxva la 1ere 405 : c *** valeur de la ligne .xxva valeur de la matrice. 406 : c *** linc(i) donne le numero de la colonne du ieme elem de xxva 407 : c 408 : c ... na = nb de lignes concernant le noeud ino ... 409 : NA = IDER-IPRE+1 410 : c ... llvnul = somme cumulée des llvva ... 411 : LLVNUL=LLVNUL+LLVVA 412 : c 413 : SEGINI,LLIGN 414 : c 415 : NBINO=INO 416 : IF (IJK.EQ.2) THEN 417 : DO 116 JJ=1,NTTMAI(/1) 418 : IF (NTTMAI(JJ).EQ.INO) GOTO 117 419 : 116 CONTINUE 420 : GOTO 118 421 : 117 CONTINUE 422 : NBINO=NNOE+JJ 423 : ENDIF 424 : 118 CONTINUE 425 : ILIGN(NBINO)=LLIGN 426 : NBA=0 427 : c 428 : c ... boucle sur les ddl du noeud ino ... 429 : DO 120 JPA=1,NA 430 : c ... iiin = n° de ligne du ddl n° jpa ... 431 : IIIN=IPRE+JPA -1 432 : c 433 : c ... que l'on met au bon endroit dans immmm (faisant partie de llign) ... 434 : c 435 : IMMMM(JPA)=IIIN 436 : c 437 : c 438 : c ... ippo(i)+1 = adresse du début (dans xxva et linc) des données relatives au ddl n° i ... 439 : IPPO(JPA)=NBA 440 : c 441 : c ... boucle sur les termes dans la ligne ... 442 : DO 121 IPAK = 1,MTRA(JPA) 443 : c ... iunpak = n° de la colonne du terme n° ipak ... 444 : IUNPAK=NTRA(IPAK,JPA) 445 : c ... pour les termes mis dans llign on efface l'information sur la position dans xtra ... 446 : LTRA(IUNPAK,JPA)=0 447 : c ... compteur ++ ... 448 : NBA=NBA+1 449 : c ... le n° de la colonne va dans linc ... 450 : LINC(NBA)=IUNPAK 451 : c ... transfert du xtra (segment tratra) vers xxva (segment llign) ... 452 : XXVA(NBA)=XTRA(IPAK,JPA) 453 : c ... les termes diagonaux vont dans diag (segment mdiag) ... 454 : IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA) 455 : 121 CONTINUE 456 : 120 CONTINUE 457 : c ... le pointeur vers la fin de la dernière ligne ... 458 : IPPO(NA+1)= NBA 459 : NJMAX= 0 460 : c ... boucle sur les lignes relatives au noeud ino ... 461 : DO 126 JPA=IPRE,IDER 462 : c ... on met le n° du noeud dans ipno (segment milign) ... 463 : IPNO(JPA)=NBINO 464 : c ... ipde et ipdf : début et fin des données relatives au ddl n° jpa dans xxva et linc ... 465 : IPDE=IPPO(JPA-IPRE+1)+1 466 : IPDF=IPPO(JPA-IPRE+2) 467 : c ... lpa = n° mini de la colonne avec des termes non nuls dans la ligne n° jpa ... 468 : LPA=JPA 469 : DO 155 JHT=IPDE,IPDF 470 : LPA=MIN(LPA,LINC(JHT)) 471 : 155 CONTINUE 472 : c ... on le met dans ldeb (segment llign) ... 473 : LDEB(JPA-IPRE+1)=LPA 474 : c ... nna = longueur de la "skyline" ... 475 : NNA= JPA- LPA +1 476 : c ... njmax = profil cumulé sur un noeud ... 477 : NJMAX=NJMAX+NNA 478 : 126 CONTINUE 479 : c 480 : c ... njtot = profil total ... 481 : NJTOT=NJTOT+NJMAX 482 : c ... ijmax = (profil / noeud) maxi ... 483 : IF(IJMAX.LT.NJMAX) IJMAX=NJMAX 484 : SEGDES,LLIGN 485 : IPRE=IDER+1 486 : 100 CONTINUE 487 : c 488 : 999 CONTINUE 489 : c 490 : c ... fin de la boucle sur les noeuds ... 491 : SEGSUP TRATRA 492 : c 493 : c **** on reprend toute les matrices contenant les multiplicateurs 494 : c **** pour multiplier tous leurs termes par une norme attachee 495 : c **** a chaque multiplicateur.puis on les assemble. 496 : c 497 : IENMU=0 498 : c ... attention ! ici commence une boucle cachée (exécutée avec un goto 375) ... 499 : 375 IENMU1 = IENMU 500 : IENMU=0 501 : c ... boucle sur les rigidités qui calcule le nombre de matrices de blocages ... 502 : DO 376 I=1,NNR 503 : IF(.NOT.INOMUL(I)) IENMU=IENMU+1 504 : 376 CONTINUE 505 : c 506 : c ... s'il n'y en a pas on saute cette partie du code ... 507 : IF( IENMU.EQ.0) GO TO 3750 508 : c 509 : c ... mimik contient les noms des variables primales ... 510 : MIMIK=IIMIK 511 : SEGACT,MIMIK 512 : c ... boucle sur les rigidités ... 513 : DO 11 I=1,NNR 514 : c ... si celle si n'est pas une matrice de blocage on passe à la suivante ... 515 : IF(INOMUL(I)) GO TO 11 516 : c 517 : DESCR=IRIGEL(3,I) 518 : c ... n3 = nb de variables primales ... 519 : N3=LISINC(/2) 520 : COER=COERIG(I) 521 : MELEME=IRIGEL(1,I) 522 : INCTRA=INCTRR(I) 523 : XMATRI=IRIGEL(4,I) 524 : c ... n2 = nombre d'éléments dans le support géométrique ... 525 : N2=NUM(/2) 526 : c ... À quoi correspond ce cas ? (pas de matrices élémentaires) ... 527 : IF (RE(/3).EQ.0) THEN 528 : INOMUL(I)=.TRUE. 529 : GOTO 11 530 : ENDIF 531 : * XMATRI=IMATTT(1) 532 : * SEGACT,XMATRI 533 : c ... n1 = nombre de variables duales ... 534 : c ... pourquoi va-t-on chercher n3 dans descr et n1 dans re ? ... 535 : N1=RE(/1) 536 : SEGINI,RA 537 : c ... boucle sur les éléments ... 538 : DO 14 IEL=1,N2 539 : c ... boucle sur les inconnues ... 540 : DO 15 ICO=1,N3 541 : c ... ija = n° local du noeud-support de l'inconnue n° ico ... 542 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 543 : c ... ijb = n° de l'inconnue ... 544 : IJB=INCTRA(ICO) 545 : c ... ival contient la correspondance entre le n° local du ddl 546 : c et le n° d'équation (de ligne) correspondant ... 547 : IVAL(ICO)=INCPO(IJB,IJA) 548 : 15 CONTINUE 549 : * DDDD= 1D-31 550 : DDDD= xpetit 551 : c ... dmax = max des valeurs absolues des termes diagonaux correspondants 552 : c au ddl de l'élément n° iel ... 553 : DMAX=DDDD 554 : DO 19 ICO=1,N3 555 : DMAX=MAX(DMAX,ABS(DIAG(IVAL(ICO)))) 556 : 19 CONTINUE 557 : c 558 : c ... aux fins d'eviter des problemes dans la decomposition 559 : c 560 : * IF(DMAX.LE.1D-30.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1) GO TO 377 561 : IF(DMAX.LE.xpetit.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1) GO TO 377 562 : * IF(DMAX.LE.1D-30) DMAX = 1.D0 563 : IF(DMAX.LE.xpetit) DMAX = 1.D0 564 : DMAX=DMAX*1.5D0 565 : c 566 : * XMATRI=IMATTT(IEL) 567 : * SEGACT,XMATRI 568 : c ... dmaxy = maximum des termes dans la première colonne (les 2 premiers exclus) ... 569 : DMAXY=0.D0 570 : DO 821 ICO=3,N1 571 : DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel))) 572 : 821 CONTINUE 573 : c 574 : IF( IIMPI. EQ.9022 ) WRITE(IOIMP,7398) DMAX 575 : 7398 FORMAT('FACTEUR MULTIPLICATIF DE NORME=',E12.5) 576 : c 577 : c ... copie de la matrice élémentaire fois dmax*coer dans ra ... 578 : DO 21 ICO=1,N1 579 : DO 2110 IKO=1,N1 580 : RA(ICO,IKO)=RE(ICO,IKO,iel)*DMAX*COER 581 : 2110 CONTINUE 582 : 21 CONTINUE 583 : c ... la sous-matrice de dimension 2 est multipliée par dmaxy ... 584 : * commenté car pose un problème dans excite et ne sert sans doute à rien 585 : * write (6,*) ' dmaxy dans assem5 ', dmaxy 586 : xcoef= 1D0 587 : RA(1,1)=RA(1,1)*xcoef 588 : RA(2,1)=RA(2,1)*xcoef 589 : RA(1,2)=RA(1,2)*xcoef 590 : RA(2,2)=RA(2,2)*xcoef 591 : c ... mise à dmax des dnor correspondant aux multiplicateurs ... 592 : DO 22 ICO=1,2 593 : c-old dnor(ival(ico))=dmax 594 : DNOR(IVAL(ICO))=DNOR(IVAL(ICO))*DMAX 595 : 22 CONTINUE 596 : c 597 : c ... boucle sur les variables primales ... 598 : DO 24 ICO=1,N3 599 : c ... ino = n° local du noeud ... 600 : INO=INUINV(NUM(NOELEP(ICO),IEL)) 601 : c ... io = n° de ligne du ddl n° ico ... 602 : IO=IVAL(ICO) 603 : if (ico.eq.1) io1=io 604 : if (ico.eq.2) io2=io 605 : c ... cas des noeuds non totalement maitres ... 606 : IF(IO.GT.NBNNMA) THEN 607 : DO 1161 JJ=1,NTTMAI(/1) 608 : IF (NTTMAI(JJ).EQ.INO) GOTO 1171 609 : 1161 CONTINUE 610 : GOTO 1181 611 : 1171 CONTINUE 612 : INO=NNOE+JJ 613 : 1181 CONTINUE 614 : ENDIF 615 : c ... llign contient les lignes liées aux noeud ino ... 616 : LLIGN=ILIGN(INO) 617 : SEGACT,LLIGN*MOD 618 : c ... on rajoute le terme diagonal au diag ... 619 : DIAG(IO)=DIAG(IO)+RA(ICO,ICO) 620 : c ... immmm contient les n°s des lignes des ddl ... 621 : DO 132 JLIJ=1,IMMMM(/1) 622 : JLIJ1=JLIJ 623 : c ... on est censé trouver le bon n° de ligne et aller au 133 ... 624 : IF( IMMMM(JLIJ).EQ.IO) GO TO 133 625 : 132 CONTINUE 626 : c 627 : IF(IIMPI.EQ.1524) WRITE(IOIMP,7354) 628 : 7354 FORMAT( ' PREMIERE ERREUR 5') 628 : CALL ERREUR(5) 629 : RETURN 630 : c 631 : 133 CONTINUE 632 : c 633 : c ... deuxième niveau de boucle sur les variables primales ... 634 : DO 26 IRO=1,N3 635 : c ... ia = n° de colonne du ddl n° iro ... 636 : IA=IVAL(IRO) 637 : c ... cas symétrique !!! ... 638 : IF(IA.GT.IO) GO TO 26 639 : c 640 : c ... jld et jlt sont égaux au début et à la fin (dans xxva et linc) 641 : c des données relatives au ddl n° iro ... 642 : JLT=IPPO(JLIJ1+1) 643 : JLD=IPPO(JLIJ1)+1 644 : c 645 : c ... on cherche le bon n° de colonne dans linc ... 646 : DO 134 JL=JLD,JLT 647 : JL1=JL 648 : IF(LINC(JL).EQ.IA) GO TO 135 649 : 134 CONTINUE 650 : c 651 : IF(IIMPI.EQ.1524) WRITE(IOIMP,7355) 652 : 7355 FORMAT( ' DEUXIEME ERREUR 5') 652 : CALL ERREUR(5) 653 : RETURN 654 : c 655 : 135 CONTINUE 656 : c ... on rajoute le coefficient dans ra au xxva correspondant ... 657 : XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO) 658 : 26 CONTINUE 659 : SEGDES,LLIGN 660 : 24 CONTINUE 661 : * on stocke dans ittr les couples de LX 662 : ittr(io1)=io2 663 : ittr(io2)=io1 664 : * SEGDES,XMATRI 665 : 14 CONTINUE 666 : c ... après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ... 667 : INOMUL(I)=.TRUE. 668 : 377 CONTINUE 669 : SEGSUP,RA 670 : 11 CONTINUE 671 : c ... fin de la boucle sur les rigidités ... 672 : GO TO 375 673 : c ... fin de la boucle cachée ... 674 : 3750 CONTINUE 675 : c 676 : DO 18 IK=1,NNR 677 : INCTRA=INCTRR(IK) 678 : SEGSUP,INCTRA 679 : 18 CONTINUE 680 : * pv on desactive tout 681 : NNR=IRIGEL(/2) 682 : DO 2 IRI=1,NNR 683 : DESCR=IRIGEL(3,IRI) 684 : SEGDES DESCR 685 : IPT1=IRIGEL(1,IRI) 686 : SEGDES IPT1 687 : xMATRI=IRIGEL(4,IRI) 688 : SEGDES xMATRI 689 : 2 CONTINUE 690 : INTERR(1)=NJTOT 691 : c 692 : IF(IIMPI.GE.1) WRITE(IOIMP,4821) LLVNUL,NJTOT 693 : 4821 FORMAT('Nbre de valeurs non nulles à l assemblage ',I9,/ 694 : & 'Nbre de valeurs à l assemblage ',I9) 695 : c 696 : c 697 : SEGDES,SNTT 698 : SEGSUP,INCTRR 699 : SEGDES,MDIAG 700 : SEGDES,MIMIK 701 : SEGDES,MDNOR 702 : SEGSUP,ITOPO 703 : SEGSUP,IITOP 704 : SEGDES,MILIGN 705 : SEGSUP,INUINV 706 : SEGDES,MMATRI 707 : MMMTRI=MMATRI 708 : SEGDES,MINCPO 709 : SEGSUP,IVAL,ITRA,JNOMUL 710 : IF(NORINC.NE.0) THEN 710 : CALL ASSEM0(MRIGID,2,INWUIT) 711 : SEGDES,MRIGID 712 : ELSE 713 : SEGDES,MRIGID 714 : ENDIF 715 : RETURN 716 : END 717 : c 718 : c 719 : c 720 : c 721 : c 722 : 723 : 724 : 725 : 726 : 727 : 728 : 729 : 730 : 731 : 732 : 733 : 734 : 735 : 736 :
© Cast3M 2003 - All rights reserved.
Disclaimer