1 : C LDMT2 SOURCE JC220346 13/12/16 21:15:56 7884 2 : SUBROUTINE LDMT2(ITRAV1,ITOPO1,INUIN1,IMINI1,MMMTRI,IPO1,INCTR1 3 : & ,INCTR2,IITOP1,TRSUP,MENAGE,LDIAG,IITOPB, 4 : & ITOPOB,IPODD,njtot,INORMU) 5 : 6 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7 : C 8 : C **** SUBROUTINE POUR FAIRE L'ASSEMBLAGE DE MATRICES 9 : C 10 : C EN ENTREE: 11 : C **** ITRAV1 : POINTEUR OBJET MRIGIDITE 12 : C **** ITOPO1 : POINTEUR SEGMENT DE TRAVAIL ITOPO ( VOIR ASSEM1) 13 : C **** IITOP1 : POINTEUR SEGMENT DE TRAVAIL IITOP ( VOIR ASSEM1) 14 : C **** INUIN1 : POINTEUR SEGMENT DE TRAVAIL INUINV(VOIR ASSEM1) 15 : C **** IMINI1 : POINTEUR SEGMENT DE TRAVAIL IMINI (VOIR ASSEM1) 16 : C **** IPO1 : POINTEUR SEGMENT DE TRAVAIL IPOS (VOIR ASSEM1) 17 : C **** MMMTRI : POINTEUR OBJET MATRICE TRIANGULARISEE (NON MODIFIE) 18 : C (VOIR SMMATRI) 19 : C **** INORMU : =0 si on ne veut pas normaliser les multiplicateurs 20 : C de Lagrange, <>0 sinon 21 : C 22 : C Appelée par : LDMT1 23 : C 24 : C Auteur : Michel BULIK 25 : C 26 : C Date : Printemps '95 27 : C 28 : C Langage : ESOPE + FORTRAN77 29 : C 30 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 31 : IMPLICIT INTEGER(I-N) 32 : -INC CCREEL 33 : -INC CCOPTIO 34 : -INC SMELEME 35 : -INC SMRIGID 36 : -INC SMMATRI 37 : 38 : SEGMENT,INUINV(NNGLOB) 39 : SEGMENT,ITOPO(IENNO) 40 : SEGMENT,ITOPOB(IENNO) 41 : SEGMENT,IITOP(NNOE+1) 42 : SEGMENT,IITOPB(NNOE+1) 43 : SEGMENT,IMINI(INC) 44 : SEGMENT,IPOS(NNOE1) 45 : SEGMENT,IPOD(NNOE1) 46 : SEGMENT,INCTRR(NIRI) 47 : SEGMENT,INCTRD(NIRI) 48 : SEGMENT,INCTRA(NLIGRE) 49 : SEGMENT,INCTRB(NLIGRE) 50 : SEGMENT,IPV(NNOE) 51 : C 52 : C **** CES TABLEAUX SERVENT AU REPERAGE DE LA MATRICE POUR L'ASSEMBLAG 53 : C **** IL SERONT TOUS SUPPRIMES EN FIN D'ASSEMBLAGE. 54 : C 55 : SEGMENT,IVAL(NNN) 56 : SEGMENT,ITRA(NNN,2) 57 : SEGMENT TRATRA 58 : REAL*8 XTRA(INCRED,INCDIF) 59 : INTEGER LTRA(INC,INCDIF) 60 : INTEGER NTRA(INCRED,INCDIF) 61 : INTEGER MTRA(INCDIF) 62 : ENDSEGMENT 63 : 64 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 65 : C 66 : C **** IVAL(I)=J : LA I EME LIGNE D'UNE PETITE MATRICE S'ASSEMBLE 67 : C DANS LA J EME DE LA GRANDE. 68 : C **** ITRAV(I,1)=J : LA IEME INCONNUE DU NOEUD EN COURS D'ASSEMBLAGE 69 : C ET QUI SE TROUVE DANS LA PETITE MATRICE SE TROUVE 70 : C EN J EME POSITION DE LA PETITE MATRICE. 71 : C **** ITRAV(I,2) : LA IEME INCONNUE DU NOEUD EN COURS D'ASSEMBLAGE 72 : C PRESENT DANS LA PETITE MATRICE EST EN JEME 73 : C POSITION DANS LA GRANDE 74 : C 75 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 76 : 77 : SEGMENT,RA(N1,N1)*D 78 : SEGMENT JNOMUL 79 : LOGICAL INOMUL(NNR) 80 : ENDSEGMENT 81 : REAL*8 DMAX,COER,DDDD,DMAXY 82 : LOGICAL NOMUL 83 : LOGICAL TRSUP,TEST,MENAGE,LDIAG 84 : 85 : 86 : SAVE IPASG 87 : DATA IPASG/0/ 88 : * 89 : * en cas de normalisation des variables. 90 : * 91 : C ... On a vérifié dans RESOU que dans le cas non-symétrique soit NORINC=NORIND=0, 92 : C soit les deux sont non-nuls, mais on se méfie quand même ... 93 : inwuit=0 94 : * write (6,*) ' NORINC NORIND LDIAG',norinc,norind,ldiag 95 : IF(NORINC.NE.0.AND.NORIND.NE.0 .AND. LDIAG) then 95 : CALL ASSEM0(ITRAV1,1,inwuit) 96 : endif 97 : * 98 : * PV ON ACTIVE UNE FOIS POUR TOUTES LES MELEME DESCR... DE LA RIGIDITE 99 : * ON EN PROFITE POUR CREER INOMUL 100 : C 101 : C **** RECHERCHE DE LA DIMENSION MAX DE IVAL,ET SEGINI DE IVAL ET ITRA 102 : C 103 : INCTRR=INCTR1 104 : SEGACT,INCTRR 105 : C ... dans le cas asymétrique on se doit de faire la même chose avec les vd ... 106 : INCTRD=INCTR2 107 : SEGACT,INCTRD 108 : 109 : MRIGID=ITRAV1 110 : SEGACT,MRIGID 111 : NNR=IRIGEL(/2) 112 : C ... NNR = équivalent à NRIGEL (nb. de rigidités élémentaires) ... 113 : 114 : NNN=0 115 : SEGINI JNOMUL 116 : DO 1 IRI=1,NNR 117 : 118 : INCTRA=INCTRR(IRI) 119 : SEGACT INCTRA 120 : C ... dans le cas asymétrique on se doit de faire la même chose avec les vd ... 121 : INCTRB=INCTRD(IRI) 122 : SEGACT INCTRB 123 : 124 : DESCR=IRIGEL(3,IRI) 125 : SEGACT, DESCR 126 : 127 : IPT1=IRIGEL(1,IRI) 128 : SEGACT IPT1 129 : 130 : xMATRI=IRIGEL(4,IRI) 131 : SEGACT xMATRI 132 : 133 : C ... NA = nb de variables primales ... 134 : NA=LISINC(/2) 135 : C ... NNN = max des nb de variables primales ... 136 : NNN=MAX(NA,NNN) 137 : 138 : C ... INOMUL (NOn MULtiplicateurs ?) dit s'il ne s'agit pas des CL ... 139 : INOMUL(IRI)=.TRUE. 140 : IF(IPT1.ITYPEL.EQ.22) INOMUL(IRI)=.FALSE. 141 : 142 : 1 CONTINUE 143 : C ... IVAL a la taille NNN ... 144 : SEGINI,IVAL 145 : C ... ITRA A LA TAILLE (NNN,2) ... 146 : SEGINI,ITRA 147 : C 148 : C **** ACTIVATION DES SEGMENTS DE TRAVAILS ET DE MMATRI 149 : C 150 : IF(.NOT.LDIAG) THEN 151 : ITOPO=ITOPO1 152 : IITOP=IITOP1 153 : IPOS=IPO1 154 : ELSE 155 : ITOPO=ITOPOB 156 : IITOP=IITOPB 157 : IPOS = IPODD 158 : ENDIF 159 : SEGACT,IITOP 160 : SEGACT,ITOPO 161 : SEGACT,IPOS 162 : * write(6,*) 'ipos', ( ipos(IU),IU=1,ipos(/1)) 163 : * write(6,*) 'itopo', ( itopo(IU),IU=1,itopo(/1)) 164 : * write(6,*) 'iitop', ( iitop(IU),IU=1,iitop(/1)) 165 : INUINV=INUIN1 166 : SEGACT,INUINV 167 : SEGACT,IPOS 168 : C ... NNOE = nombre de noeuds impliqués ... 169 : NNOE=IPOS(/1)-1 170 : C ... INC = dimension de la matrice ... 171 : INC=IPOS(NNOE+1) 172 : MMATRI=MMMTRI 173 : SEGACT,MMATRI*MOD 174 : IF(LDIAG) THEN 175 : SEGINI,MDIAG 176 : IDIAG=MDIAG 177 : ELSE 178 : MDIAG=IDIAG 179 : SEGACT,MDIAG 180 : ENDIF 181 : SEGINI,MILIGN 182 : IF(TRSUP) THEN 183 : IILIGS=MILIGN 184 : ELSE 185 : IILIGN=MILIGN 186 : ENDIF 187 : MINCPO=IINCPO 188 : SEGACT,MINCPO 189 : C ... INCDIF = nb de variables primales différentes ... 190 : INCDIF=INCPO(/1) 191 : MIPO1=IDUAPO 192 : SEGACT,MIPO1 193 : casym(OK) ... on va prendre en compte aussi les vd pour déterminer 194 : c la taille de TRATRA ... 195 : INCDID=MIPO1.INCPO(/1) 196 : INCDIF=MAX(INCDIF,INCDID) 197 : 198 : MIMIK=IIMIK 199 : SEGACT,MIMIK 200 : 201 : C ... IPV est de taille NNOE ... 202 : SEGINI IPV 203 : 204 : INCRED=0 205 : DO 80 INO =1,NNOE 206 : ICOMPT=0 207 : C ... MAXELE = nb d'éléments qui contiennent le noeud No INO ... 208 : MAXELE = (IITOP(INO+1)-IITOP(INO))/2 209 : C ... Boucle sur ces éléments ... 210 : DO 81 IELE=1,MAXELE 211 : 212 : C ... IIU = pointeur dans IITOP ... 213 : IIU=IITOP(INO) + IELE + IELE -2 214 : 215 : C ... IEL = numero de l'élément dans son maillage ... 216 : IEL=ITOPO(IIU) 217 : 218 : C ... IRI = numero du maillage (dans IRIGEL) de cet élément ... 219 : IRI=ITOPO(IIU+1) 220 : 221 : C ... On va chercher ce maillage dans IRIGEL ... 222 : MELEME=IRIGEL(1,IRI) 223 : 224 : C ... Boucle sur les noeuds de cet élément ... 225 : DO 83 I=1,NUM(/1) 226 : 227 : C ... IP = numéro local du I-ème noeud de l'élément ... 228 : IP=INUINV(NUM(I,IEL)) 229 : 230 : C ... La première condition = les "connections" sont symétriques ... 231 : IF (IP.GT.INO) GOTO 83 232 : C ... Le noeud IP a déjà vu ce noeud ... 233 : IF (IPV(IP).EQ.INO) GOTO 83 234 : 235 : C ... Ces deux opérations se font si (IP <= INO et IPV(IP) != INO) ... 236 : C ... IPV = N° max des noeuds "connectés" au noeud IP ... 237 : 238 : IPV(IP)=INO 239 : C ... ICOMPT = nombre de noeuds connectés (de numéro < ou =) au noeud INO ... 240 : ICOMPT=ICOMPT+1 241 : 83 CONTINUE 242 : 81 CONTINUE 243 : C ... INCRED = maximum des ICOMPT sur tous les noeuds impliqués ... 244 : 245 : INCRED=MAX(INCRED,ICOMPT) 246 : cdbg write(*,*) 'LDMT2 : Noeud ',INO,', ICOMPT = ',ICOMPT, 247 : cdbg & ' -> INCRED = ',INCRED 248 : 80 CONTINUE 249 : 250 : cdbg segprt,ipv 251 : SEGSUP IPV 252 : 253 : INCRED=INCRED*INCDIF 254 : C ... La taille de TRATRA dépend de INC, INCRED et INCDIF ... 255 : SEGINI TRATRA 256 : C 257 : C ... Coefficients de normalisation ... 258 : C ----------------------------- 259 : 260 : LLVNUL=0 261 : C ... On n'initialise plus IJMAX ici - ceci se fait un niveau au dessus ... 262 : C-OLD IJMAX=0 263 : 264 : C ... La taille de MDNOR est INC (nombre d'équations) ... 265 : MDNOR=0 266 : IF(LDIAG) THEN 267 : SEGINI MDNOR 268 : IDNORM=MDNOR 269 : C ... On a vérifié dans RESOU que dans le cas non-symétrique soit NORINC=NORIND=0, 270 : C soit les deux sont non-nuls ... 271 : SEGINI,MDNO1 272 : IDNORD=MDNO1 273 : DO 20 I=1,INC 274 : DNOR(I)=1.D0 275 : MDNO1.DNOR(I)=1.D0 276 : 20 CONTINUE 277 : ENDIF 278 : 279 : C ... Si une normalisation est imposée, les valeurs dans MDNOR et MDNO1 seront modifiées par ASNS3 ... 280 : midua=iidua 281 : IF(NORINC.NE.0 .AND. NORIND.NE.0 .AND. LDIAG) 281 : & CALL ASNS3(MDNOR,MIMIK,MINCPO,mdno1,midua,mipo1) 282 : 283 : C **** BOUCLE *100* SUR LES NUMEROS DE NOEUDS QUE L'ON ASSEMBLE 284 : C ------------------------------------------------------------- 285 : 286 : DO 100 INO=1,NNOE 287 : cdbg write(*,*) 'LDMT2: ----------------------------------' 288 : cdbg write(*,*) 'LDMT2: Noeud N° ',INO 289 : DO 101 IIT=1,INCDIF 290 : MTRA(IIT)=0 291 : 101 CONTINUE 292 : 293 : C ... Attention à l'égalité des IPOS et IPOD qui est testée dans ASNS1 ... 294 : 295 : C ... IPRE = N° de la première ligne (ou colonne) concernée par le noeud INO ... 296 : IPRE=IPOS(INO )+1 297 : C ... IDER = N° de la dernière ligne (ou colonne) concernée par le noeud INO ... 298 : IDER=IPOS(INO+1) 299 : cdbg write(*,*) 'LDMT2: Lignes de ',IPRE,' jusqu''à ',IDER 300 : LLVVA=0 301 : C 302 : C **** BOUCLE *99* SUR LES ELEMENTS TOUCHANT LE NOEUD INO 303 : C POUR LES ELEMNTS MULTIPLICATEUR ON NE FAIT PAS 304 : C L'ASSEMBLAGE 305 : C 306 : C ... MAXELE = nb d'éléments qui contiennent le noeud No INO ... 307 : MAXELE= (IITOP(INO+1) -IITOP(INO))/2 308 : DO 99 IELE=1,MAXELE 309 : cdbg write(*,*) 'LDMT2: - - - - - - - - - - - -' 310 : cdbg write(*,*) 'LDMT2: Elément N° ',IELE 311 : IIU=IITOP(INO) + IELE + IELE - 2 312 : C ... IEL = numero de l'élément dans son maillage ... 313 : IEL=ITOPO(IIU) 314 : C ... IRI = numero du maillage (dans IRIGEL) de cet élément ... 315 : IRI=ITOPO(IIU+1) 316 : C ... MELEME = pointeur vers ce maillage ... 317 : MELEME=IRIGEL(1,IRI) 318 : C ... DESCR = pointeur vers le segment DESCR concerné ... 319 : DESCR=IRIGEL(3,IRI) 320 : C ... INCTRA contient les indices (dans IMIK et IHAR) des descriptions 321 : C des DDL correspondants ... 322 : casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre INCTRR ou INCTRD ... 323 : INCTRA=INCTRR(IRI) 324 : INCTRB=INCTRD(IRI) 325 : C ... xMATRI = tableau de pointeurs vers les matrices élémentaires ... 326 : xMATRI=IRIGEL(4,IRI) 327 : C ... COER = coefficient multiplicateur ... 328 : COER=COERIG(IRI) 329 : C ... NOn MULtiplicateur ? ... 330 : NOMUL=INOMUL(IRI) 331 : C 332 : C **** NOMUL =.FALSE. IL EXISTE UN MULTIPLICATEUR 333 : C **** INITIALISATION DE IVAL. IVAL(I)=J VEUT DIRE QUE 334 : C **** LA I EME LIGNE DE LA PETITE MATRICE S'ASSEMBLE DANS 335 : C **** LA J EME DE LA GRANDE MATRICE. 336 : C 337 : C ... Maintenant, dans la boucle 96 (nouvelle par rapport à ASSEM2) on 338 : C va parcourir les lignes (colonnes) pour trouver celles, qui 339 : C s'appuyent sur le noeud INO, cas inférieur (supérieur) ... 340 : 341 : NA=0 342 : IF(TRSUP) THEN 343 : NIN=LISINC(/2) 344 : ELSE 345 : NIN=LISDUA(/2) 346 : ENDIF 347 : DO 96 ICO=1,NIN 348 : IF(TRSUP) THEN 349 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 350 : IJB=INCTRA(ICO) 351 : ELSE 352 : IJA=INUINV(NUM(NOELED(ICO),IEL)) 353 : IJB=INCTRB(ICO) 354 : ENDIF 355 : C ... Si on a trouvé le bon noeud ... 356 : IF(IJA.EQ.INO) THEN 357 : C ... On incrémente NA (le nombre de DDL connectés au noeud INO ?) ... 358 : NA=NA+1 359 : cdbg write(*,*) 'LDMT2: IJA == INO => NA devient ',NA 360 : ITRA(NA,1)=ICO 361 : IF(TRSUP) THEN 362 : ITRA(NA,2)=INCPO(IJB,IJA) 363 : ELSE 364 : ITRA(NA,2)=MIPO1.INCPO(IJB,IJA) 365 : ENDIF 366 : ENDIF 367 : 96 CONTINUE 368 : 369 : C ... Dans la boucle 98 on va remplir le tableau IVAL avec les numéros 370 : C de colonnes (lignes) pour tous les DDLs primaux (duaux) de la matrice, 371 : C ceci concerne le triangle inférieur (supérieur) ... 372 : casym(OK) ... la boucle devra se faire soit sur les duales soit sur les primales ... 373 : 374 : IF(TRSUP) THEN 375 : NIN=LISDUA(/2) 376 : ELSE 377 : NIN=LISINC(/2) 378 : ENDIF 379 : C ... Boucle sur les variables primales (duales) de l'élément ... 380 : DO 98 ICO=1,NIN 381 : IF(TRSUP) THEN 382 : IJA=INUINV(NUM(NOELED(ICO),IEL)) 383 : IJB=INCTRB(ICO) 384 : IVAL(ICO)=MIPO1.INCPO(IJB,IJA) 385 : cdbg write(*,*) 'LDMT2: vd locale N°',ICO,'-> vd globale N°',IJB 386 : cdbg write(*,*) 'LDMT2: => N° d''équation = ',IVAL(ICO) 387 : ELSE 388 : C ... IJA = numéro local du noeud-support du DDL ... 389 : casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre NOELEP ou NOELED ... 390 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 391 : C ... IJB = indice du DDL ... 392 : casym(OK) ... INCTRA ou INCTRB ??? ... 393 : IJB=INCTRA(ICO) 394 : C ... On met dans IVAL(ICO) le N° de la colonne correspondant au noeud IJA 395 : C et DDL No IJB ... 396 : casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre MINCPO ou MIPO1 (?) ... 397 : IVAL(ICO)=INCPO(IJB,IJA) 398 : cdbg write(*,*) 'LDMT2: vp locale N°',ICO,'-> vp globale N°',IJB 399 : cdbg write(*,*) 'LDMT2: => N° de la colonne = ',IVAL(ICO) 400 : ENDIF 401 : 98 CONTINUE 402 : 403 : * XMATRI=IMATTT(IEL) 404 : * SEGACT,XMATRI 405 : cdbg segprt,xmatri 406 : C 407 : C **** BOUCLE *95* SUR LES INCONNUES DE LA PETITE MATRICE 408 : C 409 : DO 95 INCC=1,NA 410 : C ... INCO = le N° de la ligne (cas inf.) ou colonne (cas sup.) 411 : C du DDL dual (primal) N° INCC du noeud INO ... 412 : casym(OK) ... ici il faudra choisir en fonction de TRSUP quoi prendre ... 413 : INCO=ITRA(INCC,2) 414 : 415 : C ... Boucle sur les DDL's du noeud INO ... 416 : DO 90 IK=1,NIN 417 : C ... IO = N° de la colonne (la ligne) correspondant au DDL N° IK de l'élément 418 : casym(OK) ... ici il faudra choisir en fonction de TRSUP quoi prendre ... 419 : IO=IVAL(IK) 420 : 421 : C ... Si IO > INCO on ne fait plus rien dans cette boucle ... 422 : C ... Car on ne parcourt que la moitié de la matrice !!! ... 423 : IF(IO.GT.INCO) GO TO 90 424 : 425 : C ... ILOC = N° de ligne relatif / IPRE du noeud INO ... 426 : ILOC=INCO-IPRE+1 427 : 428 : C ... JJ = N° (dans l'élément) de la ligne (de la colonne) correspondant au 429 : C DDL N° INCC du noeud INO ... 430 : casym(OK) ... attention à ce qui se passe ici ... !!! 431 : JJ=ITRA(INCC,1) 432 : 433 : C ... Ci-dessous on stockera les quelques lignes de la rigidité concernant le 434 : C noeud INO dans le segment TRATRA. Les coefficients vont dans XTRA, le deuxième 435 : C indice (ILOC) est le numéro relatif de la ligne, le premier (IMTT) n'a rien à voir avec le 436 : C numéro de la colonne. Celui-ci est stocké dans NTRA(IMTT,ILOC). L'information inverse 437 : C (IMTT en fonction de IO) se trouve dans LTRA(IO,ILOC). On met dans MTRA le nombre 438 : C de termes différents assemblés par ligne ... 439 : C ... ILTT doit contenir l'indice IMMTT en fonction du N° de la colonne; s'il est 440 : C nul, le terme n'a pas encore été assemblé, il faudra donc initialiser ... 441 : ILTT= LTRA(IO,ILOC) 442 : IF(ILTT.EQ.0) THEN 443 : C ... LLVVA est initialisé à 0 au début de la boucle 100 et sert à compter les 444 : C termes assemblés ... 445 : LLVVA=LLVVA+1 446 : C ... Les MTRA sont mis à 0 au début de la boucle 100 ... 447 : C ... Ici on incrémente MTRA(ILOC) et on met la nouvelle valeur dans IMMTT, 448 : C qui sera le premier indice dans différents tableaux de TRATRA ... 449 : IMMTT=MTRA(ILOC)+1 450 : MTRA(ILOC)=IMMTT 451 : C ... puis on initialise le nouveau terme ... 452 : XTRA(IMMTT,ILOC)=0.D0 453 : C ... son N° de colonne (cas inf.) ou ligne (cas supérieur) ... 454 : NTRA(IMMTT,ILOC)=IO 455 : C ... et l'information inverse ... 456 : LTRA(IO,ILOC)=IMMTT 457 : ILTT=IMMTT 458 : ENDIF 459 : C ... Si ce n'est pas une condition aux limites, on assemble ... 460 : IF(NOMUL) THEN 461 : C ... Attention ! Dans XTRA 1er indice est lié à la colonne, deuxième est le numéro 462 : C de la ligne, tandis que dans RE c'est l'inverse ... 463 : IF(TRSUP) THEN 464 : * write(6,*) ' iltt iloc ' , iltt,iloc 465 : * write(6,*) 'ik jj iel re coer ',ik,jj,iel,re(ik,jj,iel),coer 466 : XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IK,JJ,iel)*COER 467 : cdbg write(*,*) 'LDMT2: On rajoute le terme L',IK,', C',JJ,' soit', 468 : cdbg & RE(IK,JJ),' à la col',IPRE+ILOC-1,' et lig', 469 : cdbg & NTRA(ILTT,ILOC) 470 : ELSE 471 : * write(6,*) ' iltt iloc ' , iltt,iloc 472 : * write(6,*) 'ik jj iel re coer ',ik,jj,iel,re(jj,ik,iel),coer 473 : XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(JJ,IK,iel)*COER 474 : cdbg write(*,*) 'LDMT2: On rajoute le terme L',JJ,', C',IK,' soit', 475 : cdbg & RE(JJ,IK),' à la lig',IPRE+ILOC-1,' et col', 476 : cdbg & NTRA(ILTT,ILOC) 477 : ENDIF 478 : 479 : ENDIF 480 : 90 CONTINUE 481 : 95 CONTINUE 482 : * SEGDES,XMATRI 483 : 99 CONTINUE 484 : C 485 : C *** COMPACTAGE DES LIGNES, EN MEME TEMPS CALCUL DE IJMAX QUI SERA 486 : C *** LA DIMENSION MAX D'UN SEGMENT LIGN. 487 : C *** LE SEGMENT ASSOCIE A UNE LIGNE (SEGMENT LLIGN) EST DE LA FORME : 488 : C *** IMMMM(NA) PERMET DE SAVOIR SI UN MOUVEMENT D'ENSEMBLE SUR LA 489 : C *** LIGNE EXISTE. IPPO(NA+1) DONNE LA POSITION DANS XXVA LA 1ERE 490 : C *** VALEUR DE LA LIGNE .XXVA VALEUR DE LA MATRICE. 491 : C *** LINC(I) DONNE LE NUMERO DE LA COLONNE DU IEME ELEM DE XXVA 492 : C 493 : C ... NA = nb de lignes concernant le noeud INO ... 494 : NA = IDER-IPRE+1 495 : C ... LLVNUL = somme cumulée des LLVVA ... 496 : LLVNUL=LLVNUL+LLVVA 497 : SEGINI,LLIGN 498 : ILIGN(INO)=LLIGN 499 : NBA=0 500 : * write(6,*) ' ider ipre llvnul ino ',ider,ipre,llvnul,ino 501 : C ... Boucle sur les DDL du noeud INO ... 502 : DO 120 JPA=1,NA 503 : C ... IIIN = N° de ligne du DDL N° JPA ... 504 : IIIN=IPRE+JPA -1 505 : C ... que l'on met au bon endroit dans IMMMM (faisant partie de LLIGN) ... 506 : IMMMM(JPA)=IIIN 507 : C ... IPPO(I)+1 = adresse du début (dans XXVA et LINC) des données relatives au DDL N° I ... 508 : IPPO(JPA)=NBA 509 : C ... Boucle sur les termes dans la ligne ... 510 : DO 121 IPAK = 1,MTRA(JPA) 511 : C ... IUNPAK = N° de la colonne du terme N° IPAK ... 512 : IUNPAK=NTRA(IPAK,JPA) 513 : C ... Pour les termes mis dans LLIGN on efface l'information sur la position dans XTRA ... 514 : LTRA(IUNPAK,JPA)=0 515 : C ... Compteur ++ ... 516 : NBA=NBA+1 517 : C ... Le N° de la colonne va dans LINC ... 518 : LINC(NBA)=IUNPAK 519 : * write(6,*) ' nba ipak jpa, ', nba ,ipak,jpa 520 : C ... Transfert du XTRA (segment TRATRA) vers XXVA (segment LLIGN) ... 521 : XXVA(NBA)=XTRA(IPAK,JPA) 522 : C ... Les termes diagonaux vont dans DIAG (segment MDIAG) ... 523 : IF(LDIAG) THEN 524 : * write(6,*) ' iiin nba xxva(nba)', iiin,nba,xxva(nba) 525 : IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA) 526 : ENDIF 527 : 121 CONTINUE 528 : 120 CONTINUE 529 : C ... Le pointeur vers la fin de la dernière ligne ... 530 : IPPO(NA+1)= NBA 531 : NJMAX= 0 532 : C ... Boucle sur les lignes relatives au noeud INO ... 533 : DO 126 JPA=IPRE,IDER 534 : C ... On met le N° du noeud dans IPNO (segment MILIGN) ... 535 : IPNO(JPA)=INO 536 : C ... IPDE et IPDF : début et fin des données relatives au DDL N° JPA dans XXVA et LINC ... 537 : IPDE=IPPO(JPA-IPRE+1)+1 538 : IPDF=IPPO(JPA-IPRE+2) 539 : C ... LPA = N° mini de la colonne avec des termes non nuls dans la ligne N° JPA ... 540 : LPA=JPA 541 : DO 155 JHT=IPDE,IPDF 542 : LPA=MIN(LPA,LINC(JHT)) 543 : 155 CONTINUE 544 : C ... On le met dans LDEB (segment LLIGN) ... 545 : LDEB(JPA-IPRE+1)=LPA 546 : C ... NNA = longueur de la "skyline" ... 547 : NNA= JPA- LPA +1 548 : C ... NJMAX = profil cumulé sur un noeud ... 549 : NJMAX=NJMAX+NNA 550 : 126 CONTINUE 551 : C ... NJTOT = profil total ... 552 : NJTOT=NJTOT+NJMAX 553 : C ... IJMAX = (profil / noeud) maxi ... 554 : IF(IJMAX.LT.NJMAX) IJMAX=NJMAX 555 : SEGDES,LLIGN 556 : 100 CONTINUE 557 : C ... Fin de la boucle sur les noeuds ... 558 : SEGSUP TRATRA 559 : C 560 : C **** ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS 561 : C **** POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE 562 : C **** A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE. 563 : C 564 : IENMU=0 565 : C ... ATTENTION ! Ici commence une boucle cachée (exécutée avec un GOTO 375) ... 566 : 375 IENMU1 = IENMU 567 : IENMU=0 568 : C ... Boucle sur les rigidités qui calcule le nombre de matrices de blocages ... 569 : DO 376 I=1,NNR 570 : IF(.NOT.INOMUL(I)) IENMU=IENMU+1 571 : 376 CONTINUE 572 : 573 : C ... S'il n'y en a pas on saute cette partie du code ... 574 : IF( IENMU.EQ.0) GO TO 3750 575 : 576 : C ... MIMIK contient les noms des variables primales ... 577 : MIMIK=IIMIK 578 : SEGACT,MIMIK 579 : C ... Boucle sur les rigidités ... 580 : DO 11 I=1,NNR 581 : C ... Si celle si n'est pas une matrice de blocage on passe à la suivante ... 582 : IF(INOMUL(I)) GO TO 11 583 : 584 : DESCR=IRIGEL(3,I) 585 : C ... N3 = nb de variables primales ... 586 : N3=LISINC(/2) 587 : COER=COERIG(I) 588 : MELEME=IRIGEL(1,I) 589 : INCTRA=INCTRR(I) 590 : xMATRI=IRIGEL(4,I) 591 : C ... N2 = nombre d'éléments dans le support géométrique ... 592 : N2=NUM(/2) 593 : C ... À quoi correspond ce cas ? (pas de matrices élémentaires) ... 594 : IF (re(/3).EQ.0) THEN 595 : INOMUL(I)=.TRUE. 596 : GOTO 11 597 : ENDIF 598 : * XMATRI=IMATTT(1) 599 : * SEGACT,XMATRI 600 : C ... N1 = nombre de variables duales ... 601 : C ... Pourquoi va-t-on chercher N3 dans DESCR et N1 dans RE ? ... 602 : N1=RE(/1) 603 : SEGINI,RA 604 : C ... Boucle sur les éléments ... 605 : DO 14 IEL=1,N2 606 : 607 : C ... Boucle sur les inconnues ... 608 : DO 15 ICO=1,N3 609 : C ... IJA = N° local du noeud-support de l'inconnue N° ICO ... 610 : IJA=INUINV(NUM(NOELEP(ICO),IEL)) 611 : C ... IJB = N° de l'inconnue ... 612 : IJB=INCTRA(ICO) 613 : C ... IVAL contient la correspondance entre le N° local du DDL 614 : C et le N° d'équation (de ligne) correspondant ... 615 : IVAL(ICO)=INCPO(IJB,IJA) 616 : 15 CONTINUE 617 : C JCARDO => pour les fous qui veulent se passer de la 618 : C normalisation des mult. de Lagrange :) 619 : IF (INORMU.EQ.0) THEN 620 : DMAX=1.D0 621 : DMAXY=DMAX 622 : ELSE 623 : * DDDD= 1D-31 624 : DDDD= xpetit 625 : C ... DMAX = max des valeurs absolues des termes diagonaux correspondants 626 : C au DDL de l'élément N° IEL ... 627 : DMAX=DDDD 628 : IF((.NOT.LDIAG).AND.(IDIAG.EQ.0)) THEN 629 : * write(*,*) 'Erreur interne dans LDMT2 : MDIAG pas encore', 630 : * & ' initialisé !!!' 631 : RETURN 632 : ENDIF 633 : if (n3.lt.3) then 634 : do ico=1,n3 635 : DMAX=MAX(DMAX,ABS(DIAG(IVAL(ICO)))) 636 : enddo 637 : endif 638 : DO 19 ICO=3,N3 639 : DMAX=MAX(DMAX,ABS(DIAG(IVAL(ICO)))) 640 : 19 CONTINUE 641 : 642 : C ... AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION 643 : IF(IIMPI.EQ.1524) WRITE(IOIMP,7391) DMAX,DDDD,IENMU, 644 : & IENMU1,I,IEL 645 : 7391 FORMAT(' DMAX DDDD IENMU IENMU1 I IEL',2E12.5,4I3) 646 : 647 : * IF(DMAX.LE.1D-30.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1) GOTO 377 648 : IF(DMAX.LE.xpetit.AND.IENMU.NE.IENMU1 649 : & .AND.IEL.EQ.1) GOTO 377 650 : * IF(DMAX.LE.1D-30) DMAX = 1.D0 651 : IF(DMAX.LE.xpetit) DMAX = 1.D0 652 : DMAX=DMAX*1.5D0 653 : 654 : * XMATRI=IMATTT(IEL) 655 : * SEGACT,XMATRI 656 : C ... DMAXY = maximum des termes dans la première colonne (les 2 premiers exclus) ... 657 : DMAXY=0.D0 658 : if (n1.lt.3) then 659 : do ico=1,n1 660 : DMAXY=MAX(DMAXY,ABS( RE(ICO,1,iel))) 661 : enddo 662 : endif 663 : DO 821 ICO=3,N1 664 : DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel))) 665 : 821 CONTINUE 666 : DMAX = DMAX / DMAXY 667 : ENDIF 668 : IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7398) DMAX 669 : 7398 FORMAT(' facteur multiplicatif de norme ',e12.5) 670 : 671 : C ... Copie de la matrice élémentaire fois DMAX*COER dans RA ... 672 : DO 21 ICO=1,N1 673 : DO 2110 IKO=1,N1 674 : RA(ICO,IKO)=RE(ICO,IKO,iel)*DMAX*COER 675 : 2110 CONTINUE 676 : 21 CONTINUE 677 : C ... La sous-matrice de dimension 2 est multipliée par DMAXY ... 678 : RA(1,1)=RA(1,1)*DMAXY 679 : RA(2,1)=RA(2,1)*DMAXY 680 : RA(1,2)=RA(1,2)*DMAXY 681 : RA(2,2)=RA(2,2)*DMAXY 682 : 683 : C ... Mise à DMAX des DNOR correspondant aux multiplicateurs ... 684 : IF(LDIAG) THEN 685 : DO 22 ICO=1,2 686 : DNOR(IVAL(ICO))=DNOR(IVAL(ICO))*DMAX 687 : mdno1.DNOR(IVAL(ICO))=mdno1.DNOR(IVAL(ICO))*DMAX 688 : 22 CONTINUE 689 : ENDIF 690 : CMB ... Je suppose que l'on s'en servira pour multiplier les composantes 691 : C du second membre pour que les déplacements imposés soient corrects ... 692 : 693 : C ... Boucle sur les variables primales ... 694 : DO 24 ICO=1,N3 695 : C ... INO = N° local du noeud ... 696 : INO=INUINV(NUM(NOELEP(ICO),IEL)) 697 : C ... IO = N° de ligne du DDL N° ICO ... 698 : IO=IVAL(ICO) 699 : if (ico.eq.1) io1=io 700 : if (ico.eq.2) io2=io 701 : C ... LLIGN contient les lignes liées aux noeud INO ... 702 : LLIGN=ILIGN(INO) 703 : SEGACT,LLIGN*MOD 704 : C ... On rajoute le terme diagonal au DIAG ... 705 : IF(LDIAG) DIAG(IO)=DIAG(IO)+RA(ICO,ICO) 706 : C ... IMMMM contient les N°s des lignes des DDL ... 707 : DO 132 JLIJ=1,IMMMM(/1) 708 : JLIJ1=JLIJ 709 : C ... On est censé trouver le bon N° de ligne et aller au 133 ... 710 : IF( IMMMM(JLIJ).EQ.IO) GO TO 133 711 : 132 CONTINUE 712 : 713 : IF(IIMPI.EQ.1524) WRITE(IOIMP,7354) 714 : 7354 FORMAT( ' PREMIERE ERREUR 5') 714 : CALL ERREUR(5) 715 : RETURN 716 : 717 : 133 CONTINUE 718 : 719 : C ... Deuxième niveau de boucle sur les variables primales ... 720 : DO 26 IRO=1,N3 721 : C ... IA = N° de colonne du DDL N° IRO ... 722 : IA=IVAL(IRO) 723 : C ... Cas symétrique !!! ... 724 : IF(IA.GT.IO) GO TO 26 725 : 726 : C ... JLD et JLT sont égaux au début et à la fin (dans XXVA et LINC) 727 : C des données relatives au DDL N° IRO ... 728 : JLT=IPPO(JLIJ1+1) 729 : JLD=IPPO(JLIJ1)+1 730 : 731 : C ... On cherche le bon N° de colonne dans LINC ... 732 : DO 134 JL=JLD,JLT 733 : JL1=JL 734 : IF(LINC(JL).EQ.IA) GO TO 135 735 : 134 CONTINUE 736 : 737 : IF(IIMPI.EQ.1524) WRITE(IOIMP,7355) 738 : 7355 FORMAT( ' DEUXIEME ERREUR 5') 738 : CALL ERREUR(5) 739 : RETURN 740 : 741 : 135 CONTINUE 742 : C ... On rajoute le coefficient dans RA au XXVA correspondant ... 743 : XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO) 744 : 26 CONTINUE 745 : SEGDES,LLIGN 746 : 24 CONTINUE 747 : * on stocke dans ittr les couples de LX 748 : ittr(io1)=io2 749 : ittr(io2)=io1 750 : * SEGDES,XMATRI 751 : 14 CONTINUE 752 : C ... Après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ... 753 : INOMUL(I)=.TRUE. 754 : 377 CONTINUE 755 : SEGSUP,RA 756 : 11 CONTINUE 757 : C ... Fin de la boucle sur les rigidités ... 758 : GO TO 375 759 : C ... Fin de la boucle cachée ... 760 : 3750 CONTINUE 761 : 762 : C ... C'est l'ménaage finaale ... 763 : DO 18 IK=1,NNR 764 : INCTRA=INCTRR(IK) 765 : IF(MENAGE) THEN 766 : SEGSUP,INCTRA 767 : ELSE 768 : SEGDES,INCTRA 769 : ENDIF 770 : 18 CONTINUE 771 : * PV ON DESACTIVE TOUT 772 : NNR=IRIGEL(/2) 773 : DO 2 IRI=1,NNR 774 : DESCR=IRIGEL(3,IRI) 775 : SEGDES DESCR 776 : IPT1=IRIGEL(1,IRI) 777 : SEGDES IPT1 778 : xMATRI=IRIGEL(4,IRI) 779 : SEGDES xMATRI 780 : 2 CONTINUE 781 : INTERR(1)=NJTOT 781 : c-old IF(IPASG.EQ.0) CALL ERREUR(-278) 781 : IF(MENAGE.AND.IPASG.EQ.0) CALL ERREUR(-278) 782 : IPASG=1 783 : 784 : IF(IIMPI.EQ.1457) WRITE(IOIMP,4821) LLVNUL,NJTOT 785 : 4821 FORMAT(' NB DE VALEURS NON NULLES DANS LA MATRICE ',I9,/ 786 : & ' NB DE VALEURS DANS LA MATRICE ',I9) 787 : SEGSUP,IITOP,ITOPO 788 : IF(MENAGE) THEN 789 : C SEGSUP,IITOP 790 : ccc SEGSUP,IMINI 791 : SEGSUP,INUINV 792 : C SEGSUP,ITOPO 793 : SEGSUP,INCTRR,INCTRD 794 : ELSE 795 : C SEGDES,IITOP 796 : ccc SEGDES,IMINI 797 : SEGDES,INUINV 798 : C SEGDES,ITOPO 799 : SEGDES,INCTRR 800 : ENDIF 801 : * write(6,*) ' diag' 802 : * write(6,*) ( diag(iou),iou=1,diag(/1)) 803 : SEGDES,MDIAG 804 : SEGDES,MIMIK 805 : IF (MDNOR.NE.0) SEGDES,MDNOR 806 : SEGDES,MILIGN 807 : SEGDES,MMATRI 808 : MMMTRI=MMATRI 809 : SEGDES,MINCPO 810 : SEGDES,MIPO1 811 : SEGDES,IPOS 812 : SEGSUP,IVAL,ITRA,JNOMUL 813 : IF(NORINC.NE.0 .AND. MENAGE) THEN 813 : CALL ASSEM0(MRIGID,2,inwuit) 814 : SEGDES,MRIGID 815 : ELSE 816 : SEGDES,MRIGID 817 : ENDIF 818 : RETURN 819 : END 820 : 821 : 822 : 823 : 824 : 825 : 826 : 827 : 828 : 829 : 830 : 831 : 832 : 833 : 834 : 835 : 836 : 837 : 838 : 839 :
© Cast3M 2003 - All rights reserved.
Disclaimer