1 : C RAYN SOURCE PV 13/04/12 21:15:51 7756 2 : SUBROUTINE RAYN 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : 6 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7 : C C 8 : C But : Calculer unbe matrice de conductivité due aux échanges C 9 : C par rayonnement C 10 : C C 11 : C Syntaxe : CND1 = RAYN MO1 CHTEM CHRAY (CONSSB) ; C 12 : C C 13 : C MO1 : modèle C 14 : C CHTEM : temperature moyenne (MCHAML) C 15 : C CHRAY : matrice de rayonnement (MCHAML) C 16 : C CONSSB : constante de Stefan (REEL) C 17 : C C 18 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 19 : 20 : -INC CCOPTIO 21 : -INC CCREEL 22 : *- 23 : -INC SMCOORD 24 : -INC SMELEME 25 : -INC SMCHAML 26 : -INC SMRIGID 27 : -INC SMMODEL 28 : -INC SMINTE 29 : 30 : POINTEUR CHAMTM.MCHELM,CHAMR.MCHELM 31 : POINTEUR MELBMA.MELEME,MELCOQ.MELEME 32 : POINTEUR MAPBMA.MELEME,MAPCOQ.MELEME 33 : POINTEUR MAPCOM.MELEME 34 : 35 : SEGMENT ,INFOEL 36 : LOGICAL KCOQ(N1),KQUAD(N1) 37 : ENDSEGMENT 38 : 39 : SEGMENT,GL2LOC 40 : INTEGER IEQUIV(3,IECART) 41 : ENDSEGMENT 42 : 43 : SEGMENT/INFO/(INFELL(JG)),INFO1.INFO,INFO2.INFO 44 : SEGMENT LESINF 45 : INTEGER LINFO(NBSZEL) 46 : ENDSEGMENT 47 : 48 : SEGMENT LESAIJ 49 : REAL*8 AIJ(NBNNMX,NBELMX,NBSZEL) 50 : ENDSEGMENT 51 : 52 : DIMENSION T3(2) 53 : C ... Attention !!! 20 correspond ici au nombre maxi de noeuds par 54 : C élément 55 : C ... XEL - coordonnées d'origine, XCH - coordonnées choisies pour 56 : C VPAST (ou VPAST2), 57 : C XTR - coordonnées transformées par VCORL1 ... 58 : DIMENSION XEL(3,20),XCH(3,3),XTR(3,20) 59 : C ... BPSS - matrice de transformation ... 60 : DIMENSION BPSS(3,3) 61 : DIMENSION SHP(6,20) 62 : 63 : C-------- Fin des déclarations ------------------------------- 64 : 65 : CHARACTER*8 TYPE 66 : DATA INTTYP /3/ 67 : DATA CONSSB /5.67D-8/ 68 : 69 : C----------- Fin des DATA ------------------------------------ 70 : 71 : C ... Lecture ... 72 : 73 : TYPE='MMODEL ' 73 : CALL LIROBJ(TYPE,IRET,1,IRETOU) 74 : IF(IRETOU.EQ.0) THEN 75 : MOTERR(1:8)=TYPE 75 : CALL ERREUR(37) 76 : RETURN 77 : ENDIF 78 : IF(IRETOU.EQ.1) MMODEL=IRET 79 : SEGACT MMODEL 80 : 81 : TYPE='MCHAML ' 81 : CALL LIROBJ(TYPE,IRET,1,IRETOU) 82 : IF(IRETOU.EQ.0) THEN 83 : MOTERR(1:8)=TYPE 83 : CALL ERREUR(37) 84 : RETURN 85 : ENDIF 86 : IF(IRETOU.EQ.1) MCHEL1=IRET 87 : 88 : TYPE='MCHAML ' 88 : CALL LIROBJ(TYPE,IRET,1,IRETOU) 89 : IF(IRETOU.EQ.0) THEN 90 : MOTERR(1:8)=TYPE 90 : CALL ERREUR(37) 91 : RETURN 92 : ENDIF 93 : IF(IRETOU.EQ.1) MCHEL2=IRET 94 : 95 : SEGACT, MCHEL1, MCHEL2 96 : IF((MCHEL1.TITCHE).EQ.'MATRICE DE RAYONNEMENT') THEN 97 : CHAMR = MCHEL1 98 : CHAMTM= MCHEL2 99 : ELSEIF((MCHEL2.TITCHE).EQ.'MATRICE DE RAYONNEMENT') THEN 100 : CHAMR = MCHEL2 101 : CHAMTM= MCHEL1 102 : ELSE 102 : CALL ERREUR(21) 103 : RETURN 104 : ENDIF 105 : SEGDES, MCHEL1, MCHEL2 106 : 106 : CALL LIRREE(XVAL,0,IRETOU) 107 : IF(IRETOU.EQ.1) CONSSB=XVAL 108 : 109 : C ... Verification que les supports du MMODEL et des 110 : C MCHAML sont bien les mêmes ... 111 : 111 : CALL RAYN1(MMODEL,CHAMR) 111 : CALL RAYN1(MMODEL,CHAMTM) 112 : 113 : C CHAMTM = pointeur vers le MCHAML de température moyenne 114 : C CHAMR = pointeur vers le MCHAML de matrices de rayonnement 115 : C MMODEL = pointeur vers le modèle 116 : 117 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 118 : C ... Recherche du support géométrique de la rigidité ... C 119 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 120 : 121 : SEGACT,MMODEL 122 : 123 : C ... On créé un MELEME contenant tous les maillages-supports 124 : C du CHAMTM ... 125 : NBSZEL=KMODEL(/1) 126 : N1=NBSZEL 127 : SEGINI,INFOEL 128 : DO 1000 I=1,NBSZEL 129 : IMODEL=KMODEL(I) 130 : SEGACT,IMODEL 131 : IF(FORMOD(1).NE.'THERMIQUE'.or.matmod(2).NE.'RAYONNEMENT')THEN 131 : CALL ERREUR(21) 132 : RETURN 133 : ENDIF 134 : MELEME=IMAMOD 135 : SEGACT,MELEME 136 : 137 : IF (IDIM.EQ.3) THEN 138 : C 139 : IF ((NEFMOD.EQ.4).OR.(NEFMOD.EQ.8)) THEN 140 : C TRI3 ou QUA4 141 : KQUAD(I)=.FALSE. 142 : KCOQ(I) =.FALSE. 143 : ELSEIF ((NEFMOD.EQ.27).OR.(NEFMOD.EQ.49)) THEN 144 : C COQ3 ou COQ4 145 : KQUAD(I)=.FALSE. 146 : KCOQ(I) =.TRUE. 147 : ELSEIF ((NEFMOD.EQ.6).OR.(NEFMOD.EQ.10)) THEN 148 : C TRI6 ou QUA8 149 : KQUAD(I)=.TRUE. 150 : KCOQ(I) =.FALSE. 151 : ELSEIF ((NEFMOD.EQ.56).OR.(NEFMOD.EQ.41)) THEN 152 : C COQ6 ou COQ8 153 : KQUAD(I)=.TRUE. 154 : KCOQ(I) =.TRUE. 155 : ELSE 155 : CALL ERREUR(16) 156 : RETURN 157 : ENDIF 158 : C 159 : ELSEIF (IDIM.EQ.2) THEN 160 : C 161 : IF (NEFMOD.EQ.2) THEN 162 : C SEG2 163 : KQUAD(I)=.FALSE. 164 : KCOQ(I) =.FALSE. 165 : ELSEIF (NEFMOD.EQ.44) THEN 166 : C COQ2 167 : KQUAD(I)=.FALSE. 168 : KCOQ(I) =.TRUE. 169 : ELSEIF (NEFMOD.EQ.3) THEN 170 : C SEG3 171 : KQUAD(I)=.TRUE. 172 : KCOQ(I) =.FALSE. 173 : ELSE 173 : CALL ERREUR(16) 174 : RETURN 175 : ENDIF 176 : C 177 : ELSE 177 : CALL ERREUR(14) 178 : RETURN 179 : ENDIF 180 : 181 : C SEGDES,IMODEL 182 : 1000 CONTINUE 183 : 184 : NBMCOQ=0 185 : DO 1010 I=1,NBSZEL 186 : IF(KCOQ(I)) NBMCOQ=NBMCOQ+1 187 : 1010 CONTINUE 188 : NBMBMA=NBSZEL-NBMCOQ 189 : 190 : NBNN=0 191 : NBELEM=0 192 : NBREF=0 193 : NBSOUS=NBMBMA 194 : IF(NBMBMA.GT.0) THEN 195 : SEGINI MELBMA 196 : J=0 197 : DO 1020 I=1,NBSZEL 198 : IF(.NOT.KCOQ(I)) THEN 199 : J=J+1 200 : IMODEL=KMODEL(I) 201 : MELBMA.LISOUS(J)=IMAMOD 202 : ENDIF 203 : 1020 CONTINUE 204 : ELSE 205 : MELBMA=0 206 : ENDIF 207 : IF(IIMPI.GE.4) write(*,*) ' MELBMA',MELBMA 208 : 209 : NBNN=0 210 : NBELEM=0 211 : NBREF=0 212 : NBSOUS=NBMCOQ 213 : IF(NBMCOQ.GT.0) THEN 214 : SEGINI MELCOQ 215 : J=0 216 : DO 1030 I=1,NBSZEL 217 : IF(KCOQ(I)) THEN 218 : J=J+1 219 : IMODEL=KMODEL(I) 220 : MELCOQ.LISOUS(J)=IMAMOD 221 : ENDIF 222 : 1030 CONTINUE 223 : ELSE 224 : MELCOQ=0 225 : ENDIF 226 : IF(IIMPI.GE.4) write(*,*) ' MELCOQ',MELCOQ 227 : 228 : C ... On les transforme en maillages composés de POI1 ... 229 : IF(NBMBMA.GT.0) THEN 230 : MAPBMA=MELBMA 230 : CALL CHANGE(MAPBMA,1) 231 : 231 : c CALL ECROBJ('MAILLAGE',MAPBMA) 231 : c CALL NBNO 231 : c CALL LIRENT(NBNBMA,1,IRETOU) 232 : c IF(IERR.NE.0) GOTO 9999 233 : NBNBMA=MAPBMA.NUM(/2) 234 : ELSE 235 : MAPBMA=0 236 : NBNBMA=0 237 : ENDIF 238 : 239 : IF(NBMCOQ.GT.0) THEN 240 : MAPCOQ=MELCOQ 240 : CALL CHANGE(MAPCOQ,1) 241 : 241 : c CALL ECROBJ('MAILLAGE',MAPCOQ) 241 : c CALL NBNO 241 : c CALL LIRENT(NBNCOQ,1,IRETOU) 242 : c IF(IERR.NE.0) GOTO 9999 243 : NBNCOQ=MAPCOQ.NUM(/2) 244 : ELSE 245 : MAPCOQ=0 246 : NBNCOQ=0 247 : ENDIF 248 : 249 : C ... On détruit les maillages chapeaux ... 250 : SEGSUP,MELBMA,MELCOQ 251 : 252 : C ... On va vérifier si les deux ensembles de POI1 sont distincts ... 253 : C ... En calculant le nombre de noeuds de la différence symétrique 254 : C des deux maillages ... 255 : 256 : IF(NBNBMA*NBNCOQ.NE.0) THEN 256 : CALL ECROBJ('MAILLAGE',MAPBMA) 256 : CALL ECROBJ('MAILLAGE',MAPCOQ) 256 : CALL PRDIFF 256 : CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) 257 : IF(IIMPI.GE.4) write(*,*) ' IPT1 ',IPT1 258 : IF(IERR.NE.0) GOTO 9999 259 : 260 : C!! segprt,IPT1 261 : 261 : c CALL ECROBJ('MAILLAGE',IPT1) 261 : c CALL NBNO 261 : c CALL LIRENT(NBNDS,1,IRETOU) 262 : c IF(IERR.NE.0) GOTO 9999 263 : SEGACT IPT1 264 : NBNDS=IPT1.NUM(/2) 265 : IF(IIMPI.GE.4) THEN 266 : write(*,*) ' NBNDS ',NBNDS 267 : write(*,*) 'NBNBMA+NBNCOQ ',NBNBMA+NBNCOQ 268 : ENDIF 269 : 270 : SEGSUP,IPT1 271 : 272 : C ... Et en le comparant avec la somme des nombres de noeuds ... 273 : IF(NBNBMA+NBNCOQ.NE.NBNDS) THEN 274 : 274 : CALL ECROBJ('MAILLAGE',MAPBMA) 274 : CALL ECROBJ('MAILLAGE',MAPCOQ) 274 : CALL INTERS 274 : CALL LIROBJ('MAILLAGE',MAPCOM,1,IRETOU) 275 : IF(IERR.NE.0) GOTO 9999 276 : 276 : c CALL ECROBJ('MAILLAGE',MAPCOM) 276 : c CALL NBNO 276 : c CALL LIRENT(NBNCOM,1,IRETOU) 277 : c IF(IERR.NE.0) GOTO 9999 278 : 278 : CALL CHANGE(MAPCOM,1) 279 : NBNCOM=MAPCOM.NUM(/2) 280 : 281 : C ... Maintenant on va enlever les noeuds du MAPCOM des deux autres 282 : C maillages ... 283 : 283 : CALL ECROBJ('MAILLAGE',MAPBMA) 283 : CALL ECROBJ('MAILLAGE',MAPCOM) 283 : CALL PRDIFF 283 : CALL LIROBJ('MAILLAGE',MAPBMA,1,IRETOU) 284 : IF(IERR.NE.0) GOTO 9999 285 : NBNBMA=NBNBMA-NBNCOM 286 : 286 : CALL ECROBJ('MAILLAGE',MAPCOQ) 286 : CALL ECROBJ('MAILLAGE',MAPCOM) 286 : CALL PRDIFF 286 : CALL LIROBJ('MAILLAGE',MAPCOQ,1,IRETOU) 287 : IF(IERR.NE.0) GOTO 9999 288 : NBNCOQ=NBNCOQ-NBNCOM 289 : ELSE 290 : NBNCOM=0 291 : MAPCOM=0 292 : ENDIF 293 : 294 : ELSE 295 : NBNCOM=0 296 : MAPCOM=0 297 : ENDIF 298 : 299 : IF (IIMPI.GE.4) THEN 300 : write(*,*) 'NBNBMA = ',NBNBMA 301 : write(*,*) 'NBNCOQ = ',NBNCOQ 302 : write(*,*) 'NBNCOM = ',NBNCOM 303 : ENDIF 304 : 305 : C ... On créé le maillage - support de la rigidité ... 306 : C ... En même temps on va chercher les N° mini et maxi des noeuds ... 307 : NBNN=NBNBMA+NBNCOQ+NBNCOM 308 : NBELEM=1 309 : NBSOUS=0 310 : NBREF=NBSZEL 311 : SEGINI,MELEME 312 : ITYPEL=28 313 : 314 : IF(NBNBMA.GT.0) THEN 315 : SEGACT,MAPBMA 316 : NOMINI=MAPBMA.NUM(1,1) 317 : NOMAXI=MAPBMA.NUM(1,1) 318 : DO 1100 I=1,NBNBMA 319 : NUM(I,1)=MAPBMA.NUM(1,I) 320 : IF(MAPBMA.NUM(1,I).LT.NOMINI) NOMINI=MAPBMA.NUM(1,I) 321 : IF(MAPBMA.NUM(1,I).GT.NOMAXI) NOMAXI=MAPBMA.NUM(1,I) 322 : 1100 CONTINUE 323 : SEGDES,MAPBMA 324 : ENDIF 325 : 326 : IF(NBNCOQ.GT.0) THEN 327 : SEGACT,MAPCOQ 328 : IF(NBNBMA.EQ.0) THEN 329 : NOMINI=MAPCOQ.NUM(1,1) 330 : NOMAXI=MAPCOQ.NUM(1,1) 331 : ENDIF 332 : DO 1110 I=1,NBNCOQ 333 : NUM(I+NBNBMA,1)=MAPCOQ.NUM(1,I) 334 : IF(MAPCOQ.NUM(1,I).LT.NOMINI) NOMINI=MAPCOQ.NUM(1,I) 335 : IF(MAPCOQ.NUM(1,I).GT.NOMAXI) NOMAXI=MAPCOQ.NUM(1,I) 336 : 1110 CONTINUE 337 : SEGDES,MAPCOQ 338 : ENDIF 339 : 340 : IF(NBNCOM.GT.0) THEN 341 : SEGACT,MAPCOM 342 : IF(NBNBMA+NBNCOQ.EQ.0) THEN 343 : NOMINI=MAPCOM.NUM(1,1) 344 : NOMAXI=MAPCOM.NUM(1,1) 345 : ENDIF 346 : DO 1120 I=1,NBNCOM 347 : NUM(I+NBNBMA+NBNCOQ,1)=MAPCOM.NUM(1,I) 348 : IF(MAPCOM.NUM(1,I).LT.NOMINI) NOMINI=MAPCOM.NUM(1,I) 349 : IF(MAPCOM.NUM(1,I).GT.NOMAXI) NOMAXI=MAPCOM.NUM(1,I) 350 : 1120 CONTINUE 351 : SEGDES,MAPCOM 352 : ENDIF 353 : 354 : DO 1200 I=1,NBREF 355 : IMODEL=KMODEL(I) 356 : LISREF(I)=IMAMOD 357 : c-inutile SEGDES,IMODEL 358 : C ... On profite de cette boucle pour activer les maillages ... 359 : IPT1=IMAMOD 360 : SEGACT,IPT1 361 : 1200 CONTINUE 362 : ICOLOR(1)=7 363 : 364 : cdebug segprt,meleme 365 : SEGDES,MELEME 366 : 367 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 368 : C ... Préparation de la structure du MRIGID ... C 369 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 370 : 371 : NRIGE=7 372 : NRIGEL=1 373 : SEGINI,MRIGID 374 : MTYMAT='RIGIDITE' 375 : COERIG(1)=1.D0 376 : IRIGEL(1,1)=MELEME 377 : C ... Les noeuds des massifs ont un seul DDL : T 378 : C coques ont deux DDL : TINF TSUP 379 : C communs ont trois DDL : T TINF et TSUP 380 : NLIGRP=NBNBMA+2*NBNCOQ+3*NBNCOM 381 : NLIGRD=NLIGRP 382 : SEGINI,DESCR 383 : 384 : DO 2000 I=1,NBNBMA 385 : LISINC(I)='T ' 386 : LISDUA(I)='Q ' 387 : NOELEP(I)=I 388 : NOELED(I)=I 389 : 2000 CONTINUE 390 : 391 : DO 2010 I=1,NBNCOQ 392 : LISINC(NBNBMA+2*I-1)='TSUP' 393 : LISINC(NBNBMA+2*I )='TINF' 394 : LISDUA(NBNBMA+2*I-1)='QSUP' 395 : LISDUA(NBNBMA+2*I )='QINF' 396 : NOELEP(NBNBMA+2*I-1)=NBNBMA+I 397 : NOELEP(NBNBMA+2*I )=NBNBMA+I 398 : NOELED(NBNBMA+2*I-1)=NBNBMA+I 399 : NOELED(NBNBMA+2*I )=NBNBMA+I 400 : 2010 CONTINUE 401 : 402 : DO 2020 I=1,NBNCOM 403 : LISINC(2*NBNCOQ+NBNBMA+3*I-2)='T' 404 : LISINC(2*NBNCOQ+NBNBMA+3*I-1)='TSUP' 405 : LISINC(2*NBNCOQ+NBNBMA+3*I )='TINF' 406 : LISDUA(2*NBNCOQ+NBNBMA+3*I-2)='Q' 407 : LISDUA(2*NBNCOQ+NBNBMA+3*I-1)='QSUP' 408 : LISDUA(2*NBNCOQ+NBNBMA+3*I )='QINF' 409 : NOELEP(2*NBNCOQ+NBNBMA+3*I-2)=2*NBNCOQ+NBNBMA+I 410 : NOELEP(2*NBNCOQ+NBNBMA+3*I-1)=2*NBNCOQ+NBNBMA+I 411 : NOELEP(2*NBNCOQ+NBNBMA+3*I )=2*NBNCOQ+NBNBMA+I 412 : NOELED(2*NBNCOQ+NBNBMA+3*I-2)=2*NBNCOQ+NBNBMA+I 413 : NOELED(2*NBNCOQ+NBNBMA+3*I-1)=2*NBNCOQ+NBNBMA+I 414 : NOELED(2*NBNCOQ+NBNBMA+3*I )=2*NBNCOQ+NBNBMA+I 415 : 2020 CONTINUE 416 : 417 : IRIGEL(3,1)=DESCR 418 : 419 : cdebug segprt,descr 420 : SEGDES,DESCR 421 : 422 : NELRIG=1 423 : * SEGINI,IMATRI 424 : SEGINI,XMATRI 425 : C ... On remplira XMATRI plus tard ... 426 : * IMATTT(1)=XMATRI 427 : * SEGDES,IMATRI 428 : 429 : IRIGEL(4,1)=xMATRI 430 : IRIGEL(5,1)=NIFOUR 431 : IRIGEL(6,1)=0 432 : C ... On est dans le cas asymétrique ... 433 : IRIGEL(7,1)=2 434 : IFORIG=IFOMOD 435 : SEGDES,MRIGID 436 : 437 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 438 : C ... Calcul proprement dit de la matrice de rigidité ... C 439 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 440 : 441 : C ... On commence par préparer le segment GL2LOC ... 442 : IECART=NOMAXI-NOMINI+1 443 : SEGINI,GL2LOC 444 : 445 : C ... et le remplir ... 446 : IF(NBNBMA.GT.0) THEN 447 : SEGACT,MAPBMA 448 : DO 3000 I=1,NBNBMA 449 : IEQUIV(1,MAPBMA.NUM(1,I)-NOMINI+1)=I 450 : 3000 CONTINUE 451 : SEGDES,MAPBMA 452 : ENDIF 453 : 454 : IF(NBNCOQ.GT.0) THEN 455 : SEGACT,MAPCOQ 456 : DO 3010 I=1,NBNCOQ 457 : IEQUIV(2,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I-1 458 : IEQUIV(3,MAPCOQ.NUM(1,I)-NOMINI+1)=NBNBMA+2*I 459 : 3010 CONTINUE 460 : SEGDES,MAPCOQ 461 : ENDIF 462 : 463 : IF(NBNCOM.GT.0) THEN 464 : SEGACT,MAPCOM 465 : DO 3020 I=1,NBNCOQ 466 : IEQUIV(1,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-2 467 : IEQUIV(2,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I-1 468 : IEQUIV(3,MAPCOM.NUM(1,I)-NOMINI+1)=NBNBMA+2*NBNCOQ+3*I 469 : 3020 CONTINUE 470 : SEGDES,MAPCOM 471 : ENDIF 472 : 473 : cdebug segprt,gl2loc 474 : 475 : C ... On détruit les MAPOIN ? dont on n'aura plus besoin ? vraiment ?... 476 : *** IF(NBNBMA.GT.0) SEGSUP,MAPBMA 477 : IF(NBNCOQ.GT.0) SEGSUP,MAPCOQ 478 : IF(NBNCOM.GT.0) SEGSUP,MAPCOM 479 : 480 : C ... On commence par extraire un segment INFO pour chaque zone 481 : C élémentaire. On profite de cette boucle pour trouver le 482 : C maximum des NBNN sur toutes les zones élémentaires ... 483 : 484 : SEGINI,LESINF 485 : DO 3900 I=1,NBSZEL 486 : 487 : IMODEL=KMODEL(I) 488 : MELE=NEFMOD 489 : segact imodel 490 : if( 2+inttyp.gt.infmod(/1) ) then 490 : CALL ELQUOI(MELE,0,INTTYP,ikk,IMODEL) 491 : info = ikk 492 : else 493 : jg=16 494 : segini info 495 : do iou=1,16 496 : infell(iou)=infele(iou) 497 : enddo 498 : infell(11)=infmod(2+inttyp) 499 : endif 500 : LINFO(I)=info 501 : MELEME=IMAMOD 502 : NBNN=NUM(/1) 503 : NBEL=NUM(/2) 504 : IF(I.EQ.1) THEN 505 : NBNNMX = NBNN 506 : NBELMX = NBEL 507 : ENDIF 508 : IF(NBNN.GT.NBNNMX) NBNNMX = NBNN 509 : IF(NBEL.GT.NBELMX) NBELMX = NBEL 510 : 511 : 3900 CONTINUE 512 : 513 : C ... NBNNMX et NBELMX servent à dimensionner le segment LESAIJ 514 : C qui contiendra les termes suivants : 515 : C 516 : C ik / 517 : C a = | N sur l'élément N° k de la zone élémentaire N° i, 518 : C j / j j étant l'indice du noeud de l'élément (et de 519 : C la fonction de forme associée) 520 : C 521 : 522 : SEGINI,LESAIJ 523 : IDIM1=IDIM-1 524 : DO 3910 I=1,NBSZEL 525 : 526 : IMODEL=KMODEL(I) 527 : MELEME=IMAMOD 528 : NBNN=NUM(/1) 529 : NBEL=NUM(/2) 530 : 531 : INFO=LINFO(I) 532 : MINTE=INFELL(11) 533 : SEGACT,MINTE 534 : cdebug write(*,*) '-----------------------------' 535 : cdebug write(*,*) 'ZONE N° ',I 536 : cdebug segprt,minte 537 : NBNGAU=POIGAU(/1) 538 : 539 : DO 3915 IEL=1,NBEL 540 : 540 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL) 541 : 542 : C ... SEG2 ... 543 : IF(ITYPEL.EQ.2) THEN 544 : XCH(1,1)=XEL(1,1) 545 : XCH(2,1)=XEL(2,1) 546 : XCH(1,2)=XEL(1,2) 547 : XCH(2,2)=XEL(2,2) 548 : C ... SEG3 ... 549 : ELSE IF(ITYPEL.EQ.3) THEN 550 : XCH(1,1)=XEL(1,1) 551 : XCH(2,1)=XEL(2,1) 552 : XCH(1,2)=XEL(1,3) 553 : XCH(2,2)=XEL(2,3) 554 : C ... TRI3 ou QUA4 ... 555 : ELSE IF(ITYPEL.EQ.4 .OR. ITYPEL.EQ.8) THEN 556 : XCH(1,1)=XEL(1,1) 557 : XCH(2,1)=XEL(2,1) 558 : XCH(3,1)=XEL(3,1) 559 : XCH(1,2)=XEL(1,2) 560 : XCH(2,2)=XEL(2,2) 561 : XCH(3,2)=XEL(3,2) 562 : XCH(1,3)=XEL(1,3) 563 : XCH(2,3)=XEL(2,3) 564 : XCH(3,3)=XEL(3,3) 565 : C ... TRI6 ou QUA8 ... 566 : ELSE IF(ITYPEL.EQ.6 .OR. ITYPEL.EQ.10) THEN 567 : XCH(1,1)=XEL(1,1) 568 : XCH(2,1)=XEL(2,1) 569 : XCH(3,1)=XEL(3,1) 570 : XCH(1,2)=XEL(1,3) 571 : XCH(2,2)=XEL(2,3) 572 : XCH(3,2)=XEL(3,3) 573 : XCH(1,3)=XEL(1,5) 574 : XCH(2,3)=XEL(2,5) 575 : XCH(3,3)=XEL(3,5) 576 : ELSE 576 : CALL ERREUR(16) 577 : RETURN 578 : ENDIF 579 : 580 : IF(IDIM.EQ.2)THEN 580 : CALL VPAST2(XCH,BPSS) 581 : ELSE IF(IDIM.EQ.3) THEN 581 : CALL VPAST(XCH,BPSS) 582 : ENDIF 583 : 583 : CALL VCORL1(XEL,XTR,BPSS,NBNN) 584 : 585 : DO 3920 IPG=1,NBNGAU 586 : 587 : DO 3930 INN=1,NBNN 588 : SHP(1,INN)=SHPTOT(1,INN,IPG) 589 : SHP(2,INN)=SHPTOT(2,INN,IPG) 590 : SHP(3,INN)=SHPTOT(3,INN,IPG) 591 : 3930 CONTINUE 592 : 592 : CALL JACOBI(XTR,SHP,IDIM1,NBNN,DJAC) 593 : cdebug write(*,*) 'Zone ',I,', Element ',IEL,', PG ',IPG, 594 : cdebug & '-> J = ',DJAC 595 : 596 : IF(IFOMOD.NE.0) THEN 597 : 598 : DO 3935 INN=1,NBNN 599 : AIJ(INN,IEL,I) = AIJ(INN,IEL,I) + 600 : & SHPTOT(1,INN,IPG)*POIGAU(IPG)*DJAC 601 : 3935 CONTINUE 602 : 603 : ELSE 604 : C ... axisymetrique : calcul du rayon au point d'integration 605 : C 605 : CALL DISTRR(XEL,SHP,NBNN,RR) 606 : RR=RR*2.D0*XPI 607 : 608 : DO 3936 INN=1,NBNN 609 : AIJ(INN,IEL,I) = AIJ(INN,IEL,I) + 610 : & SHPTOT(1,INN,IPG)*POIGAU(IPG)*RR*DJAC 611 : 3936 CONTINUE 612 : ENDIF 613 : 614 : 3920 CONTINUE 615 : 616 : 3915 CONTINUE 617 : 618 : 3910 CONTINUE 619 : 620 : cdebug segprt,lesaij 621 : 622 : C ... Double boucle sur tous les éléments concernés ... 623 : C ... Tout ce qui se termine par 1 concerne la boucle extérieure : 624 : C I1, IMODE1, IPT1, INFO1, MINTE1 ... 625 : 626 : SEGACT,CHAMTM 626 : CALL TCHAMR(CHAMR,1) 627 : 628 : BILAN = 0.D0 629 : DO 4000 I1=1,NBSZEL 630 : 631 : IMODE1=KMODEL(I1) 632 : IPT1=IMODE1.IMAMOD 633 : NBNN1=IPT1.NUM(/1) 634 : 635 : C ... On va chercher la température moyenne de l'élément IEL1 ... 636 : C ... On commence par trouver les bons numéros de composantes, 637 : C puis on active les MELVAL associés ... 638 : 639 : SEGACT,CHAMTM 640 : MCHAM1=CHAMTM.ICHAML(I1) 641 : C write(*,*) ' I1 MCHAM1 = ',I1,MCHAM1 642 : SEGACT,MCHAM1 643 : NC=MCHAM1.IELVAL(/1) 644 : IF(KCOQ(I1)) THEN 645 : NCTINF=0 646 : NCTSUP=0 647 : DO 4005 I=1,NC 648 : IF(MCHAM1.NOMCHE(I).EQ.'TINF ') NCTINF=I 649 : IF(MCHAM1.NOMCHE(I).EQ.'TSUP ') NCTSUP=I 650 : 4005 CONTINUE 651 : 652 : IF(NCTINF.EQ.0) GOTO 9999 653 : MELVA1=MCHAM1.IELVAL(NCTINF) 654 : SEGACT,MELVA1 655 : 656 : IF(NCTSUP.EQ.0) GOTO 9999 657 : MELVA2=MCHAM1.IELVAL(NCTSUP) 658 : SEGACT,MELVA2 659 : ELSE 660 : NCT=0 661 : DO 4006 I=1,NC 662 : IF(MCHAM1.NOMCHE(I).EQ.'T ') NCT=I 663 : 4006 CONTINUE 664 : 665 : IF(NCT.EQ.0) GOTO 9999 666 : MELVA1=MCHAM1.IELVAL(NCT) 667 : SEGACT,MELVA1 668 : ENDIF 669 : SEGDES,MCHAM1 670 : SEGDES,CHAMTM 671 : 672 : DO 4010 IEL1=1,IPT1.NUM(/2) 673 : 674 : C ... On trouve la (les) température(s) correspondants à 675 : C l'élément IEL1 ... 676 : 677 : IF(KCOQ(I1)) THEN 678 : IF(MELVA1.VELCHE(/2).EQ.1) THEN 679 : T=MELVA1.VELCHE(1,1) 680 : ELSE 681 : T=MELVA1.VELCHE(1,IEL1) 682 : ENDIF 683 : T3(2)=T*T*T 684 : 685 : IF(MELVA2.VELCHE(/2).EQ.1) THEN 686 : T=MELVA2.VELCHE(1,1) 687 : ELSE 688 : T=MELVA2.VELCHE(1,IEL1) 689 : ENDIF 690 : T3(1)=T*T*T 691 : ELSE 692 : IF(MELVA1.VELCHE(/2).EQ.1) THEN 693 : T=MELVA1.VELCHE(1,1) 694 : ELSE 695 : T=MELVA1.VELCHE(1,IEL1) 696 : ENDIF 697 : T3(1)=T*T*T 698 : ENDIF 699 : 700 : C ... Deuxième niveau de la boucle sur les éléments ... 701 : 702 : DO 4020 I2=1,NBSZEL 703 : 704 : IMODE2=KMODEL(I2) 705 : IPT2=IMODE2.IMAMOD 706 : NBNN2=IPT2.NUM(/1) 707 : 708 : DO 4030 IEL2=1,IPT2.NUM(/2) 709 : 710 : C ... On va chercher les surfaces des éléments IEL1 et IEL2 ... 711 : 712 : MCHAML=CHAMR.ICHAML(I2) 713 : SEGACT,MCHAML 714 : NC2=IELVAL(/1) 715 : C ... Ici on suppose que la surface est toujours dernière ... 716 : MELVAL=IELVAL(NC2) 717 : SURF2=VELCHE(1,IEL2) 718 : 719 : MCHAML=CHAMR.ICHAML(I1) 720 : SEGACT,MCHAML 721 : NC1=IELVAL(/1) 722 : C ... Ici on suppose que la surface est toujours dernière ... 723 : MELVAL=IELVAL(NC1) 724 : SURF1=VELCHE(1,IEL1) 725 : 726 : C ... On a supposé que dans le champ de facteurs de forme les 727 : C composantes sont toujours dans l'ordre suivant : 728 : C MIDL SURF dans le cas des solides 729 : C SUPE INFE SURF dans le cas des coques 730 : C 731 : DO 4040 IC1=1,NC1-1 732 : MELVAL=IELVAL(IC1) 733 : MCHELM=IELCHE(1,IEL1) 734 : MCHAML=ICHAML(I2) 735 : DO 4050 IC2=1,NC2-1 736 : MELVAL=IELVAL(IC2) 737 : COEFF=VELCHE(1,IEL2) 738 : C write(*,*)' iel1 iel2 coeff ',IEL1,IEL2,COEFF 739 : 740 : DO 4060 INN1=1,NBNN1 741 : NGL1=IPT1.NUM(INN1,IEL1) 742 : NDDL1=IC1 743 : IF(KCOQ(I1)) NDDL1=NDDL1+1 744 : NLOC1 = IEQUIV(NDDL1,NGL1-NOMINI+1) 745 : DO 4070 INN2=1,NBNN2 746 : 747 : NGL2=IPT2.NUM(INN2,IEL2) 748 : NDDL2=IC2 749 : IF(KCOQ(I2)) NDDL2=NDDL2+1 750 : NLOC2 = IEQUIV(NDDL2,NGL2-NOMINI+1) 751 : 752 : C write(*,*)' in1 in2 coeff ',INN1,INN2,COEFF 753 : TERME = COEFF*AIJ(INN1,IEL1,I1) 754 : & *AIJ(INN2,IEL2,I2)*T3(IC1)*CONSSB/(SURF1*SURF2) 755 : RE(NLOC2,NLOC1,1)=RE(NLOC2,NLOC1,1) + TERME 756 : BILAN = BILAN + TERME 757 : 758 : cdebug write(*,*) 'On rajoute ',TERME,' au terme ',NLOC1,NLOC2 759 : 760 : 4070 CONTINUE 761 : 762 : 4060 CONTINUE 763 : 764 : 4050 CONTINUE 765 : 766 : 4040 CONTINUE 767 : 768 : 4030 CONTINUE 769 : 770 : 4020 CONTINUE 771 : 772 : 4010 CONTINUE 773 : SEGDES,MELVA1 774 : IF(KCOQ(I1)) SEGDES,MELVA2 775 : 776 : 4000 CONTINUE 777 : IF (IIMPI.GE.4) THEN 778 : write(*,*) 'Bilan = ',BILAN 779 : ENDIF 780 : 781 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 782 : C ... Sortie de la rigidité calculée ... C 783 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 784 : 785 : C!! segprt,xmatri 786 : SEGDES,XMATRI 786 : CALL ECROBJ('RIGIDITE',MRIGID) 787 : 788 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 789 : C ... Le ménage à la fin ... C 790 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 791 : 792 : C ... On désactive le modèle ... 793 : DO 5000 I=1,NBSZEL 794 : IMODEL=KMODEL(I) 795 : MELEME=IMAMOD 796 : SEGDES,MELEME 797 : SEGDES,IMODEL 798 : INFO=LINFO(I) 799 : SEGSUP,INFO 800 : 5000 CONTINUE 801 : SEGDES,MMODEL 802 : SEGSUP,LESINF 803 : 804 : C ... desactivation de la matrice de rayonnement 804 : CALL TCHAMR(CHAMR,0) 805 : 806 : C ... On supprime les segments de travail ... 807 : SEGSUP,INFOEL 808 : SEGSUP,GL2LOC 809 : SEGSUP,LESAIJ 810 : 811 : RETURN 812 : 813 : 9999 CONTINUE 813 : CALL ERREUR(26) 814 : RETURN 815 : END 816 : 817 : 818 : 819 : 820 :
© Cast3M 2003 - All rights reserved.
Disclaimer