1 : C MASSE4 SOURCE AM 11/06/09 21:15:01 7004 2 : SUBROUTINE MASSE4(IPMAIL,LW,LRE,IVAMAT,NMATT,IVACAR,NCARR, 3 : &IVECT,NBPGAU,IPMINT,NDDL,MELE,MFR,IPMATR,ILUMP, 4 : &ISOUS,IIPDPG,IMOD) 5 : *---------------------------------------------------------------------* 6 : * ________________________________ * 7 : * | | * 8 : * | calcul de la matrice de masse | * 9 : * |________________________________| * 10 : * * 11 : * raccords liquide/massifs,raccords liquide/coque,barre,homogenise * 12 : * cerce,joints 2-3d * 13 : * * 14 : *---------------------------------------------------------------------* 15 : * * 16 : * entrees : * 17 : * ________ * 18 : * * 19 : * ipmail pointeur sur un segment meleme * 20 : * lw dimension du tableau de travail de l'element * 21 : * lre nombre de ddl dans la matrice de masse * 22 : * ivamat pointeur sur un segment mptval pour le materiau * 23 : * nmatt nombre de composante de materiau (imat=1) * 24 : * ivacar pointeur sur un segment mptval pour les caracteri- * 25 : * stiques * 26 : * ncarr nombre de caracteristiques geometriques * 27 : * ivect flag indiquant si on a entree les axes locaux * 28 : * nbpgau nombre de point d'integration pour la masse * 29 : * ipmint pointeur sur un segment minte * 30 : * ipmin1 pointeur sur un segment minte (aux noeuds) * 31 : * nddl nombre de degre de liberte /noeud * 32 : * mele numero de l'element fini * 33 : * nbpgmi nombre de noeuds /element * 34 : * ilump =1 si l'op�ateur LUMP est appel� * 35 : * * 36 : * sorties : * 37 : * ________ * 38 : * * 39 : * ipmatr pointeur sur la matrice de masse de la sous-zone * 40 : * * 41 : *---------------------------------------------------------------------* 42 : IMPLICIT INTEGER(I-N) 43 : IMPLICIT REAL*8(A-H,O-Z) 44 : -INC CCOPTIO 45 : -INC CCHAMP 46 : -INC CCREEL 47 : -INC SMRIGID 48 : -INC SMCHAML 49 : -INC SMELEME 50 : -INC SMCOORD 51 : -INC SMINTE 52 : -INC SMMODEL 53 : -INC SMLMOTS 54 : c 55 : SEGMENT WRK1 56 : REAL*8 REL(LRE,LRE),XE(3,NBBB) 57 : ENDSEGMENT 58 : c 59 : SEGMENT WRK2 60 : REAL*8 SHPWRK(6,NBNO),BGENE(NDDL,LRE) 61 : ENDSEGMENT 62 : c 63 : SEGMENT WRK3 64 : REAL*8 WORK(LW) 65 : ENDSEGMENT 66 : * 67 : SEGMENT WRK4 68 : REAL*8 BPSS(3,3),XEL(3,NBBB) 69 : ENDSEGMENT 70 : * 71 : SEGMENT WRK5 72 : REAL*8 XGENE(NCOM,LRN) 73 : ENDSEGMENT 74 : * 75 : SEGMENT MVELCH 76 : REAL*8 VALMAT(NV1) 77 : ENDSEGMENT 78 : c 79 : SEGMENT MPTVAL 80 : INTEGER IPOS(NS),NSOF(NS) 81 : INTEGER IVAL(NCOSOU) 82 : CHARACTER*16 TYVAL(NCOSOU) 83 : ENDSEGMENT 84 : 85 : CHARACTER*8 CMATE 86 : CHARACTER*4 lesinc(7),lesdua(7) 87 : DATA lesinc/'UX','UY','UZ','RX','RY','RZ','UR'/ 88 : DATA lesdua/'FX','FY','FZ','MX','MY','MZ','FR'/ 89 : * 90 : MELEME=IPMAIL 91 : NBNN=NUM(/1) 92 : NBELEM=NUM(/2) 93 : * 94 : NV1=NMATT 95 : SEGINI,MVELCH 96 : * 97 : xMATRI=IPMATR 98 : * lval = (lre*(lre+1))/2 99 : NLIGRP=LRE 100 : NLIGRD=LRE 101 : * 102 : * introduction du point autour duquel se fait le mouvement 103 : * de la section en defo plane generalisee 104 : * 105 : IF (IFOUR.EQ.-3.AND.MFR.NE.35) THEN 106 : IP=IIPDPG 107 : SEGACT MCOORD 108 : IREF=(IP-1)*(IDIM+1) 109 : XDPGE=XCOOR(IREF+1) 110 : YDPGE=XCOOR(IREF+2) 111 : ELSE 112 : XDPGE=0.D0 113 : YDPGE=0.D0 114 : ENDIF 115 : * 116 : NHRM=NIFOUR 117 : * 118 : MINTE=IPMINT 119 : * 120 : IMODEL = IMOD 121 : CMATE = CMATEE 122 : jmat = 0 123 : iinc=0 124 : idua=0 125 : DO imat = 1 , matmod(/2) 126 : if (matmod(imat).eq.'IMPEDANCE') then 127 : jmat = imat 128 : goto 45 129 : endif 130 : ENDDO 131 : IF (mfr.eq.28) THEN 132 : jgn = 4 133 : if (ifour.eq.2) then 134 : jgm = 6 135 : segini mlmots 136 : iinc = mlmots 137 : do igm = 1,jgm 138 : mots(igm) = lesinc(igm) 139 : enddo 140 : segini mlmots 141 : idua = mlmots 142 : do igm= 1,jgm 143 : mots(igm) = lesdua(igm) 144 : enddo 145 : else if (ifour.lt.0) then 146 : jgm = 4 147 : segini mlmots 148 : iinc = mlmots 149 : mots(1) = lesinc(1) 150 : mots(2) = lesinc(2) 151 : mots(3) = lesinc(4) 152 : mots(4) = lesinc(5) 153 : segini mlmots 154 : idua = mlmots 155 : mots(1) = lesdua(1) 156 : mots(2) = lesdua(2) 157 : mots(3) = lesdua(4) 158 : mots(4) = lesdua(5) 159 : else if (ifour.eq.0) then 160 : jgm = 3 161 : segini mlmots 162 : iinc = mlmots 163 : mots(1) = lesinc(7) 164 : mots(2) = lesinc(3) 165 : mots(3) = lesinc(6) 166 : segini mlmots 167 : idua = mlmots 168 : mots(1) = lesdua(7) 169 : mots(2) = lesdua(3) 170 : mots(3) = lesdua(6) 171 : else if (ifour.eq.1) then 172 : * a faire 173 : endif 174 : ENDIF 175 : 176 : c 177 : c numero des etiquettes : 178 : c etiquettes de 1 a 98 pour traitement specifique a l element 179 : c dans la zone specifique a chaque element commencant par : 180 : c 5 continue 181 : c element 5 etiquettes 1005 2005 3005 4005 ... 182 : c 44 continue 183 : c element 44 etiquettes 1044 2044 3044 4044 ... 184 : c_______________________________________________________________________ 185 : c 186 : * CABL SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 QUA9 187 : GOTO ( 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 188 : * RAC2 RAC3 CUB8 CU20 PRI6 PR15 LIA3 LIA4 LIA6 LIA8 MULT 189 : & , 12, 99, 99, 99, 99, 99, 18, 18, 99, 99, 99 190 : * TET4 TE10 PYR5 PY13 COQ3 DKT POUT LISP FAC3 FAC4 FAC6 191 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 192 : * FAC8 LTR3 LQU4 LCU8 LPR6 LTE4 LPY5 COQ8 TUYA TUFI COQ2 193 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 194 : * POI1 BARR RACO LSU2 COQ4 LISM COF3 RES2 LSU3 LSU4 LICO 195 : & , 45, 46, 47, 99, 99, 99, 99, 99, 99, 99, 55 196 : * COQ6 CVS2 CVS3 CVT3 CVT6 CVQ4 CVQ8 THP5 TH13 THP6 TH15 197 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 198 : * THC8 TH20 ICT3 ICQ4 ICT6 ICQ8 ICC8 ICT4 ICP6 IC20 IC10 199 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 200 : * IC15 TRIP QUAP CUBP TETP PRIP TIMO JOI2 JOI3 JOT3 JOI4 201 : & , 99, 99, 99, 99, 99, 99, 99, 85, 99, 87, 88 202 : * JOI6 JOI8 LISC TRIH DST LIC4 CERC TUYO LSE2 LITU HYT3 203 : & , 99, 99, 99, 92, 99, 94, 46, 99, 99, 12, 99 204 : * HYQ4 HYT4 HYP6 HYC8 TRIS QUAS POIS FOR3 JOP3 JOP6 JOP8 205 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 206 : * POL3 POL4 POL5 POL6 POL7 POL8 POL9 PO10 PO11 PO12 PO13 207 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 208 : * PO14 BAR3 BAEX LIA2 QUAH CUBH ROT3 SEF2 TRF3 QUF4 CUF8 209 : & , 99, 46, 124, 99, 126, 127, 99, 99, 99, 99, 99 210 : * PRF6 TEF4 PYF5 MSE3 MTR6 MQU9 MC27 MP18 MT10 MP14 SEF3 211 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 212 : * TRF7 QUF9 CF27 PF21 TF15 PF19 SEG6 TR21 QU36 C216 P126 213 : & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 214 : * TE56 PY91 TRH6 215 : & , 99, 99, 157),MELE 216 : * 217 : GOTO(168,169,170,171,172),MELE-167 218 : * 219 : * JOI1 220 : GOTO(265),MELE-264 221 : 222 : c 223 : 224 : GOTO 99 225 : c_______________________________________________________________________ 226 : c 227 : c secteur de calcul pour les elements de raccord rac2 rac3 litu 228 : c liquide massif lineaire cas bidimensionnel 229 : c_______________________________________________________________________ 230 : c 231 : 12 CONTINUE 232 : IF (ILUMP .EQ. 1) GOTO 99 233 : NBBB=NBNN 234 : IF (MELE.NE.98) LW=IDIM 235 : SEGINI WRK1,WRK3 236 : DO 3012 IB=1,NBELEM 237 : c 238 : c on cherche les coordonnees de l element ib 239 : c 239 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 239 : CALL ZERO(REL,LRE,LRE) 240 : c 241 : c calcul des coefficients de normalisation 242 : c 243 : MPTVAL=IVAMAT 244 : c 245 : IF (MELE.NE.98) THEN 246 : c 247 : DO 5012 IM=1,NMATT 248 : MELVAL=IVAL(IM) 249 : IBMN=MIN(IB,VELCHE(/2)) 250 : VALMAT(IM)=VELCHE(1,IBMN) 251 : 5012 CONTINUE 252 : RHOREF=VALMAT(1) 253 : RLCAR = VALMAT(2) 254 : c 255 : CFPI= RHOREF*RLCAR 256 : c 257 : ELSE 258 : c 259 : c cas de l'element litu 260 : c 261 : DO 7012 IM=1,NMATT 262 : MELVAL=IVAL(IM) 263 : IBMN=MIN(IB,VELCHE(/2)) 264 : WORK(IM+9)=VELCHE(1,IBMN) 265 : 7012 CONTINUE 266 : ENDIF 267 : c 268 : c lecture des caracteristiques dans work 269 : c 270 : MPTVAL=IVACAR 271 : DO 4012 IC=1,NCARR 272 : IF (IVAL(IC).NE.0) THEN 273 : MELVAL=IVAL(IC) 274 : IBMN=MIN(IB,VELCHE(/2)) 275 : WORK(IC)=VELCHE(1,IBMN) 276 : ELSE 277 : WORK(IC)=0.D0 278 : ENDIF 279 : 4012 CONTINUE 280 : c 281 : c cas ou on a lu le mot vecteur 282 : c 283 : IF (IVECT.EQ.1) THEN 284 : IF (IVAL(NCARR).NE.0) THEN 285 : MELVAL=IVAL(NCARR) 286 : IBMN=MIN(IB,IELCHE(/2)) 287 : IP=IELCHE(1,IBMN) 288 : IREF=(IP-1)*(IDIM+1) 289 : DO 6112 IC=1,IDIM 290 : WORK(NCARR+IC-1)=XCOOR(IREF+IC) 291 : 6112 CONTINUE 292 : ELSE 293 : DO 6212 IC=1,IDIM 294 : WORK(NCARR+IC-1)=0 295 : 6212 CONTINUE 296 : ENDIF 297 : ENDIF 298 : c 299 : c 300 : c 301 : IF (MELE.EQ.98) THEN 301 : CALL COUMAS(REL,LRE,WORK,XE,KERRE) 302 : ELSE 302 : CALL RACMAS(NBPGAU,IFOUR,NIFOUR,IDIM,NBNN,XE,CFPI,WORK, 303 : 1 POIGAU,SHPTOT,REL,LRE) 304 : ENDIF 305 : c 306 : * SEGINI XMATRI 307 : * IMATTT(IB)=XMATRI 308 : c 309 : c remplissage de xmatri 310 : c 310 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 311 : * SEGDES XMATRI 312 : 3012 CONTINUE 313 : SEGDES xMATRI 314 : SEGSUP WRK1,WRK3,MVELCH 315 : GOTO 510 316 : c_______________________________________________________________________ 317 : c 318 : c secteur de calcul pour les elements de raccord lia3 lia4 319 : c liquide massif lineaire cas tridimensionnel 320 : c_______________________________________________________________________ 321 : c 322 : 18 CONTINUE 323 : IF (ILUMP .EQ. 1) GOTO 99 324 : NBBB=NBNN 325 : LW=IDIM 326 : SEGINI WRK1,WRK3 327 : DO 3018 IB=1,NBELEM 328 : c 329 : c on cherche les coordonnees de l element ib 330 : c 330 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 330 : CALL ZERO(REL,LRE,LRE) 331 : c 332 : c calcul des coefficients de normalisation 333 : c 334 : MPTVAL=IVAMAT 335 : DO 5018 IM=1,NMATT 336 : MELVAL=IVAL(IM) 337 : IBMN=MIN(IB,VELCHE(/2)) 338 : VALMAT(IM)=VELCHE(1,IBMN) 339 : 5018 CONTINUE 340 : RHOREF=VALMAT(1) 341 : RLCAR = VALMAT(2) 342 : c 343 : CFPI= RHOREF*RLCAR 344 : c 345 : MPTVAL=IVACAR 346 : DO 4018 IC=1,NCARR 347 : IF (IVAL(IC).NE.0) THEN 348 : MELVAL=IVAL(IC) 349 : IBMN=MIN(IB,VELCHE(/2)) 350 : WORK(IC)=VELCHE(1,IBMN) 351 : ELSE 352 : WORK(IC)=0.D0 353 : ENDIF 354 : 4018 CONTINUE 355 : c 355 : CALL LIAMAS(NBPGAU,IDIM,NBNN,NDDL,XE,CFPI,WORK,POIGAU, 356 : 1 SHPTOT,REL,LRE,IER246) 357 : IF(IER246.NE.0) THEN 357 : CALL ERREUR(IER246) 358 : SEGSUP xMATRI 359 : SEGSUP WRK1,WRK3,MVELCH 360 : GOTO 510 361 : ENDIF 362 : * SEGINI XMATRI 363 : * IMATTT(IB)=XMATRI 364 : c 365 : c remplissage de xmatri 366 : c 366 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 367 : * SEGDES XMATRI 368 : 3018 CONTINUE 369 : SEGDES xMATRI 370 : SEGSUP WRK1,WRK3,MVELCH 371 : GOTO 510 372 : c_______________________________________________________________________ 373 : c 374 : c element point (poi1) en defos planes generalisees 375 : c_______________________________________________________________________ 376 : c 377 : 45 CONTINUE 378 : IF (jmat.gt.0) THEN 379 : MPTVAL=IVAMAT 380 : MELVAL=IVAL(1) 381 : if (ival(/1).ge.1) then 382 : * melva1 = ival(1) 383 : melva1 = ival(2) 384 : else 385 : melva1 = 0 386 : endif 387 : jddl = LRE/NBPGAU 388 : DO IB = 1,NBELEM 389 : JDIAG = 0 390 : IBMN=MIN(IB,VELCHE(/2)) 391 : ig = 1 392 : igmn = MIN(IG,VELCHE(/1)) 393 : XMASS=VELCHE(IGMN,IBMN) 394 : XINER = XMASS 395 : if (melva1.gt.0) then 396 : XINER = melva1.VELCHE(IGMN,IBMN) 397 : endif 398 : do idl = 1,jddl 399 : JDIAG = JDIAG + 1 400 : RE(JDIAG,JDIAG,ib) = XMASS 401 : if (idl.gt.3) RE(JDIAG,JDIAG,ib) = XINER 402 : enddo 403 : do idl = jddl+1, LRE 404 : JDIAG = JDIAG + 1 405 : RE(JDIAG,JDIAG,ib) = XMASS 406 : if (idl.gt.jddl+3) RE(JDIAG,JDIAG,ib) = XINER 407 : enddo 408 : ENDDO 409 : SEGDES xMATRI 410 : GOTO 510 411 : ENDIF 412 : 413 : IF (MFR.EQ.26) THEN 414 : * MODAL (car goto 510 sous IMPEDANCE) 415 : DO IB = 1,NBELEM 416 : * SEGINI XMATRI 417 : * IMATTT(IB)=XMATRI 418 : 419 : MPTVAL=IVAMAT 420 : MELVAL=IVAL(2) 421 : IBMN=MIN(IB,VELCHE(/2)) 422 : XMASS=VELCHE(1,IBMN) 423 : RE(1,1,ib) = XMASS 424 : * SEGDES XMATRI 425 : ENDDO 426 : SEGDES xMATRI 427 : GOTO 510 428 : * 429 : ELSE IF (MFR.EQ.28) THEN 430 : * STATIQUE (car goto 510 sous IMPEDANCE) 431 : DO IB = 1,NBELEM 432 : * SEGINI XMATRI 433 : * IMATTT(IB)=XMATRI 434 : 435 : MPTVAL=IVAMAT 436 : MELVAL=IVAL(1) 437 : IBMN=MIN(IB,IELCHE(/2)) 438 : idepl=IELCHE(1,IBMN) 439 : MELVAL=IVAL(3) 440 : IBMN=MIN(IB,IELCHE(/2)) 441 : imade=IELCHE(1,IBMN) 441 : CALL XTY1(idepl,imade,iinc,idua,X1) 442 : re(1,1,ib) = x1 443 : * SEGDES XMATRI 444 : ENDDO 445 : SEGDES xMATRI 446 : GOTO 510 447 : ENDIF 448 : * 449 : IF(MELE.EQ.45.AND.IFOUR.NE.-3) GOTO 99 450 : NBBB=NBNN 451 : SEGINI WRK1,WRK3 452 : c 453 : c boucle de calcul pour les differents elements 454 : c 455 : DO 3045 IB=1,NBELEM 456 : c 457 : c on cherche les coordonnees de l element ib 458 : c 458 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 459 : c 460 : c on recherche rho et la section 461 : c 462 : MPTVAL=IVAMAT 463 : MELVAL=IVAL(1) 464 : IBMN=MIN(IB,VELCHE(/2)) 465 : RR=VELCHE(1,IBMN) 466 : MPTVAL=IVACAR 467 : MELVAL=IVAL(1) 468 : IBMN=MIN(IB,VELCHE(/2)) 469 : RR=RR*VELCHE(1,IBMN) 470 : c 471 : c on calcule la matrice de masse 472 : c 472 : CALL PO1MAS(XE,XDPGE,YDPGE,RR,LRE,REL) 473 : c 474 : * SEGINI XMATRI 475 : * IMATTT(IB)=XMATRI 475 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 476 : * SEGDES XMATRI 477 : 3045 CONTINUE 478 : SEGDES xMATRI 479 : SEGSUP WRK1,WRK3,MVELCH 480 : GO TO 510 481 : c_______________________________________________________________________ 482 : c 483 : c elements barre et cerce 484 : c_______________________________________________________________________ 485 : c 486 : 46 CONTINUE 487 : * 488 : IF(MELE.EQ.95.AND.IFOUR.NE.0.AND.IFOUR.NE.1) GOTO 99 489 : * 490 : NBBB=NBNN 491 : SEGINI WRK1,WRK3 492 : IF(MELE.EQ.123) THEN 493 : NCOM=NBNN 494 : LRN =LRE 495 : SEGINI WRK5 496 : ENDIF 497 : c 498 : c boucle de calcul pour les differents elements 499 : c 500 : DO 3046 IB=1,NBELEM 501 : c 502 : c on cherche les coordonnees de l element ib 503 : c 503 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 504 : c 505 : MPTVAL=IVAMAT 506 : MELVAL=IVAL(1) 507 : IBMN=MIN(IB,VELCHE(/2)) 508 : RR=VELCHE(1,IBMN) 509 : MPTVAL=IVACAR 510 : MELVAL=IVAL(1) 511 : IBMN=MIN(IB,VELCHE(/2)) 512 : RR=RR*VELCHE(1,IBMN) 513 : c 514 : c on calcule la matrice de masse 515 : c 515 : IF(MELE.EQ.46) CALL BARMAS(REL,LRE,RR,XE) 515 : IF(MELE.EQ.95) CALL CERMAS(REL,LRE,RR,XE) 515 : IF(MELE.EQ.123) CALL MASBA3(REL,LRE,RR,XE,XGENE,KERRE) 516 : IF(KERRE.NE.0) INTERR(1)=ISOUS 517 : IF(KERRE.NE.0) INTERR(2)=IB 517 : IF(MELE.EQ.123.AND.KERRE.EQ.1) CALL ERREUR(128) 518 : c 519 : * SEGINI XMATRI 520 : * IMATTT(IB)=XMATRI 521 : IF (ILUMP .EQ. 1) THEN 521 : CALL LUMP1(REL,LRE,RE(1,1,ib),IFOUR) 522 : ELSE 522 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 523 : ENDIF 524 : * SEGDES XMATRI 525 : 3046 CONTINUE 526 : SEGDES xMATRI 527 : SEGSUP WRK1,WRK3,MVELCH 528 : IF(MELE.EQ.123) SEGSUP WRK5 529 : GO TO 510 530 : c_______________________________________________________________________ 531 : c 532 : c JOINT UNIDIMENSIONNEL JOI1 533 : c_______________________________________________________________________ 534 : c 535 : 265 CONTINUE 536 : * 537 : NBBB=NBNN 538 : SEGINI WRK1,WRK3,WRK4 539 : * 540 : DO 3265 IB=1,NBELEM 541 : c 542 : MPTVAL=IVAMAT 543 : DO IC=1,NMATT 544 : IF(IVAL(IC).NE.0) THEN 545 : MELVAL=IVAL(IC) 546 : IBMN=MIN(IB,VELCHE(/2)) 547 : WORK(IC)=VELCHE(1,IBMN) 548 : ELSE 549 : WORK(IC)=0.D0 550 : ENDIF 551 : END DO 552 : * 552 : CALL MAPALU(NMATT,WORK,BPSS,IDIM) 553 : c 554 : c on calcule la matrice de masse localement 555 : c 555 : CALL JOIMAS(REL,LRE,WORK,NMATT,IDIM) 556 : c 557 : c on passe en repère global 558 : c 559 : IAW1=101 560 : IAW2=IAW1+LRE*LRE 561 : IAW3=IAW2+LRE*LRE 562 : IAW4=IAW3+LRE*LRE 562 : CALL JOIGLO(REL,BPSS,WORK(IAW1),WORK(IAW2), 563 : & WORK(IAW3),WORK(IAW4),LRE,IDIM) 564 : * 565 : IF (ILUMP .EQ. 1) THEN 565 : CALL LUMP1(REL,LRE,RE(1,1,ib),IFOUR) 566 : ELSE 566 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 567 : ENDIF 568 : 3265 CONTINUE 569 : SEGDES xMATRI 570 : SEGSUP WRK1,WRK3,WRK4,MVELCH 571 : GO TO 510 572 : c_______________________________________________________________________ 573 : c 574 : c element barre 3d excentre (baex) 575 : c_______________________________________________________________________ 576 : c 577 : 124 CONTINUE 578 : NBBB=NBNN 579 : NCOM=2 580 : LRN =LRE 581 : SEGINI WRK1,WRK3,WRK5 582 : c 583 : c boucle de calcul pour les differents elements 584 : c 585 : DO 3199 IB=1,NBELEM 586 : c 587 : c on recupere la section de l'element, ses excentrements et son 588 : c orientation. les caracteristiques sont rangees dans work 589 : c selon l'ordre suivant: sect excz excy vx vy vz 590 : c 591 : MPTVAL=IVACAR 592 : DO IC=1,NCARR 593 : IF(IVAL(IC).NE.0) THEN 594 : MELVAL=IVAL(IC) 595 : IBMN=MIN(IB,VELCHE(/2)) 596 : WORK(IC)=VELCHE(1,IBMN) 597 : ELSE 598 : WORK(IC)=0.D0 599 : ENDIF 600 : END DO 601 : SECT=WORK(1) 602 : c 603 : c xgene stocke la matrice de passage de l'element excentre 604 : c 604 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 604 : CALL MAPAEX(XE,NBNN,WORK,AL,XGENE,LRE,KERRE) 605 : IF(KERRE.NE.0) INTERR(1)=ISOUS 606 : IF(KERRE.NE.0) INTERR(2)=IB 606 : IF(KERRE.EQ.1) CALL ERREUR(128) 607 : c 608 : MPTVAL=IVAMAT 609 : MELVAL=IVAL(1) 610 : IBMN=MIN(IB,VELCHE(/2)) 611 : RR=VELCHE(1,IBMN)*SECT 612 : c 613 : c on calcule la matrice de masse 614 : c 614 : CALL BAMAEX(REL,LRE,RR,AL,XGENE) 615 : c 616 : * SEGINI XMATRI 617 : * IMATTT(IB)=XMATRI 618 : IF (ILUMP .EQ. 1) THEN 618 : CALL LUMP1(REL,LRE,RE(1,1,ib),IFOUR) 619 : ELSE 619 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 620 : ENDIF 621 : * SEGDES XMATRI 622 : 3199 CONTINUE 623 : SEGDES xMATRI 624 : SEGSUP WRK1,WRK3,WRK5,MVELCH 625 : GO TO 510 626 : c_______________________________________________________________________ 627 : c 628 : c secteur de calcul pour les elements de raccord 629 : c liquide coque cas bidimensionnel 630 : c_______________________________________________________________________ 631 : c 632 : 47 CONTINUE 633 : IF (ILUMP .EQ. 1) GOTO 99 634 : NBBB=NBNN 635 : LW=IDIM 636 : SEGINI WRK1,WRK3 637 : DO 3047 IB=1,NBELEM 638 : c 639 : c on cherche les coordonnees de l element ib 640 : c 640 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 640 : CALL ZERO(REL,LRE,LRE) 641 : c 642 : c calcul des coefficients de normalisation 643 : c 644 : MPTVAL=IVAMAT 645 : DO 5047 IM=1,NMATT 646 : MELVAL=IVAL(IM) 647 : IBMN=MIN(IB,VELCHE(/2)) 648 : VALMAT(IM)=VELCHE(1,IBMN) 649 : 5047 CONTINUE 650 : RHOREF=VALMAT(1) 651 : RLCAR = VALMAT(2) 652 : c 653 : CFPI= RHOREF*RLCAR 654 : c 655 : MPTVAL=IVACAR 656 : DO 4047 IC=1,NCARR 657 : IF (IVAL(IC).NE.0) THEN 658 : MELVAL=IVAL(IC) 659 : IBMN=MIN(IB,VELCHE(/2)) 660 : WORK(IC)=VELCHE(1,IBMN) 661 : ELSE 662 : WORK(IC)=0.D0 663 : ENDIF 664 : 4047 CONTINUE 665 : c 665 : CALL RACOMA(IFOUR,NIFOUR,IDIM,XE,CFPI,WORK,REL,LRE) 666 : c 667 : * SEGINI XMATRI 668 : * IMATTT(IB)=XMATRI 669 : c 670 : c remplissage de xmatri 671 : c 671 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 672 : * SEGDES XMATRI 673 : 3047 CONTINUE 674 : SEGDES xMATRI 675 : SEGSUP WRK3 676 : SEGSUP WRK1,MVELCH 677 : GOTO 510 678 : c_______________________________________________________________________ 679 : c 680 : c secteur de calcul pour les elements de raccord 681 : c liquide coque 3 noeuds - cas tridimensionnel 682 : c_______________________________________________________________________ 683 : c 684 : 55 CONTINUE 685 : IF (ILUMP .EQ. 1) GOTO 99 686 : NBBB=NBNN 687 : LW=IDIM 688 : SEGINI WRK1,WRK3 689 : DO 3055 IB=1,NBELEM 690 : c 691 : c on cherche les coordonnees de l element ib 692 : c 692 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 692 : CALL ZERO(REL,LRE,LRE) 693 : c 694 : c calcul des coefficients de normalisation 695 : c 696 : MPTVAL=IVAMAT 697 : DO 5055 IM=1,NMATT 698 : MELVAL=IVAL(IM) 699 : IBMN=MIN(IB,VELCHE(/2)) 700 : VALMAT(IM)=VELCHE(1,IBMN) 701 : 5055 CONTINUE 702 : RHOREF=VALMAT(1) 703 : RLCAR = VALMAT(2) 704 : c 705 : CFPI= RHOREF*RLCAR 706 : c 707 : MPTVAL=IVACAR 708 : DO 4055 IC=1,NCARR 709 : IF (IVAL(IC).NE.0) THEN 710 : MELVAL=IVAL(IC) 711 : IBMN=MIN(IB,VELCHE(/2)) 712 : WORK(IC)=VELCHE(1,IBMN) 713 : ELSE 714 : WORK(IC)=0.D0 715 : ENDIF 716 : 4055 CONTINUE 717 : c 717 : CALL LICOMA(NBPGAU,IDIM,NBNN,NDDL,XE,CFPI,WORK,POIGAU, 718 : 1 QSIGAU,ETAGAU,SHPTOT,REL,LRE,IER246) 719 : IF(IER246.NE.0) THEN 719 : CALL ERREUR(IER246) 720 : SEGSUP xMATRI 721 : SEGSUP WRK1,WRK3,MVELCH 722 : GOTO 510 723 : ENDIF 724 : * SEGINI XMATRI 725 : * IMATTT(IB)=XMATRI 726 : c 727 : c remplissage de xmatri 728 : c 728 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 729 : * SEGDES XMATRI 730 : 3055 CONTINUE 731 : SEGDES xMATRI 732 : SEGSUP WRK1,WRK3,MVELCH 733 : GOTO 510 734 : c_______________________________________________________________________ 735 : c 736 : c secteur de calcul pour les elements joints joi2 737 : c_______________________________________________________________________ 738 : c 739 : 85 CONTINUE 740 : IF (ILUMP .EQ. 1) GOTO 99 741 : NBNO=NBNN 742 : NBBB=NBNN 743 : SEGINI WRK1,WRK2,WRK4 744 : I195=0 745 : I259=0 746 : DO 3085 IB=1,NBELEM 747 : c 748 : c on cherche les coordonnees des noeuds de l element ib 749 : c 749 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 749 : CALL ZERO (REL,LRE,LRE) 750 : c 751 : c calcul des coordonnees locales de l'element 752 : c 752 : CALL JO2LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 753 : c 754 : c boucle sur les points de gauss 755 : c 756 : ISDJC=0 757 : DO 4085 IGAU=1,NBPGAU 757 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 758 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 759 : * 760 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 761 : IF(DJAC.EQ.0.) I259=IB 762 : DJAC=ABS(DJAC)*POIGAU(IGAU) 763 : MPTVAL=IVAMAT 764 : IF (IVAL(1).NE.0) THEN 765 : MELVAL=IVAL(1) 766 : IGMN=MIN(IGAU,VELCHE(/1)) 767 : IBMN=MIN(IB,VELCHE(/2)) 768 : VALMAT(1)=VELCHE(IGMN,IBMN) 769 : ELSE 770 : VALMAT(1)=0.D0 771 : ENDIF 772 : CCCCCCCCCCC DJAC=DJAC*VALMAT(1)/3.0D0 773 : C 774 : C IL FAUT DIVISER PAR 4, CE QUI CORRESPOND PLUS EXACTEMENT A DIVISER 775 : C LE B PAR 2... 776 : C 777 : DJAC=DJAC*VALMAT(1)/4.0D0 778 : c 779 : c cas axisymetrique : multiplication par le rayon de courbure 780 : c 781 : IF (IFOUR.EQ.0) THEN 782 : RAYON = 0.D0 783 : NUMSUP=NBNO/2 784 : DO 4185 IRAY=1,NUMSUP 785 : RAYON=RAYON + ( SHPTOT(1,IRAY,IGAU)*XE(1,IRAY) ) 786 : 4185 CONTINUE 787 : DJAC=DJAC*RAYON 788 : ENDIF 789 : c 789 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 790 : 4085 CONTINUE 791 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 792 : * SEGINI XMATRI 793 : * IMATTT(IB)=XMATRI 794 : c 795 : c remplissage de xmatri 796 : c 796 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 797 : * SEGDES XMATRI 798 : 3085 CONTINUE 799 : IF(I195.NE.0) INTERR(1)=I195 799 : IF(I195.NE.0) CALL ERREUR(195) 800 : IF(I259.NE.0) INTERR(1)=I259 800 : IF(I259.NE.0) CALL ERREUR(259) 801 : SEGDES xMATRI 802 : SEGSUP WRK1,WRK2,WRK4,MVELCH 803 : GOTO 510 804 : c_______________________________________________________________________ 805 : c 806 : c secteur de calcul pour les elements joints jot3 807 : c_______________________________________________________________________ 808 : c 809 : 87 CONTINUE 810 : IF (ILUMP .EQ. 1) GOTO 99 811 : NBNO=NBNN 812 : NBBB=NBNN 813 : SEGINI WRK1,WRK2,WRK4 814 : I195=0 815 : I259=0 816 : DO 3087 IB=1,NBELEM 817 : c 818 : c on cherche les coordonnees des noeuds de l element ib 819 : c 819 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 819 : CALL ZERO (REL,LRE,LRE) 820 : c 821 : c calcul des coordonnees locales de l'element 822 : c 822 : CALL JT3LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 823 : c 824 : c boucle sur les points de gauss 825 : c 826 : ISDJC=0 827 : DO 4087 IGAU=1,NBPGAU 827 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 828 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 829 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 830 : IF(DJAC.EQ.0.) I259=IB 831 : DJAC=ABS(DJAC)*POIGAU(IGAU) 832 : MPTVAL=IVAMAT 833 : IF (IVAL(1).NE.0) THEN 834 : MELVAL=IVAL(1) 835 : IGMN=MIN(IGAU,VELCHE(/1)) 836 : IBMN=MIN(IB,VELCHE(/2)) 837 : VALMAT(1)=VELCHE(IGMN,IBMN) 838 : ELSE 839 : VALMAT(1)=0.D0 840 : ENDIF 841 : DJAC=DJAC*VALMAT(1) 842 : CCCCCCCCCCC DJAC=DJAC/3.0D0 843 : C IL FAUT DIVISER PAR 4, CE QUI CORRESPOND PLUS EXACTEMENT A DIVISER 844 : C LE B PAR 2... 845 : DJAC=DJAC/4.0D0 845 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 846 : 4087 CONTINUE 847 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 848 : * SEGINI XMATRI 849 : * IMATTT(IB)=XMATRI 850 : c 851 : c remplissage de xmatri 852 : c 852 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 853 : * SEGDES XMATRI 854 : 3087 CONTINUE 855 : IF(I195.NE.0) INTERR(1)=I195 855 : IF(I195.NE.0) CALL ERREUR(195) 856 : IF(I259.NE.0) INTERR(1)=I259 856 : IF(I259.NE.0) CALL ERREUR(259) 857 : SEGDES xMATRI 858 : SEGSUP WRK1,WRK2,WRK4,MVELCH 859 : GOTO 510 860 : c_______________________________________________________________________ 861 : c 862 : c secteur de calcul pour les elements joints joi4 863 : c_______________________________________________________________________ 864 : c 865 : 88 CONTINUE 866 : IF (ILUMP .EQ. 1) GOTO 99 867 : NBNO=NBNN 868 : NBBB=NBNN 869 : SEGINI WRK1,WRK2,WRK4 870 : I195=0 871 : I259=0 872 : DO 3088 IB=1,NBELEM 873 : c 874 : c on cherche les coordonnees des noeuds de l element ib 875 : c 875 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 875 : CALL ZERO (REL,LRE,LRE) 876 : c 877 : c calcul des coordonnees locales de l'element 878 : c 878 : CALL JO4LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 879 : c 880 : c boucle sur les points de gauss 881 : c 882 : ISDJC=0 883 : DO 4088 IGAU=1,NBPGAU 883 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 884 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 885 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 886 : IF(DJAC.EQ.0.) I259=IB 887 : DJAC=ABS(DJAC)*POIGAU(IGAU) 888 : MPTVAL=IVAMAT 889 : IF (IVAL(1).NE.0) THEN 890 : MELVAL=IVAL(1) 891 : IGMN=MIN(IGAU,VELCHE(/1)) 892 : IBMN=MIN(IB,VELCHE(/2)) 893 : VALMAT(1)=VELCHE(IGMN,IBMN) 894 : ELSE 895 : VALMAT(1)=0.D0 896 : ENDIF 897 : DJAC=DJAC*VALMAT(1) 898 : CCCCCCCCCCCC DJAC=DJAC/3.0D0 899 : C IL FAUT DIVISER PAR 4, CE QUI CORRESPOND PLUS EXACTEMENT A DIVISER 900 : C LE B PAR 2... 901 : DJAC=DJAC/4.0D0 901 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 902 : 4088 CONTINUE 903 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 904 : * SEGINI XMATRI 905 : * IMATTT(IB)=XMATRI 906 : c 907 : c remplissage de xmatri 908 : c 908 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 909 : * SEGDES XMATRI 910 : 3088 CONTINUE 911 : IF(I195.NE.0) INTERR(1)=I195 911 : IF(I195.NE.0) CALL ERREUR(195) 912 : IF(I259.NE.0) INTERR(1)=I259 912 : IF(I259.NE.0) CALL ERREUR(259) 913 : SEGDES xMATRI 914 : SEGSUP WRK1,WRK2,WRK4,MVELCH 915 : GOTO 510 916 : c_______________________________________________________________________ 917 : c 918 : c secteur de calcul pour les elements joints jgi2 919 : c_______________________________________________________________________ 920 : c 921 : 170 CONTINUE 922 : IF (IFOUR.EQ.-3) NDDL=NDDL+1 923 : IF (ILUMP .EQ. 1) GOTO 99 924 : NBNO=NBNN 925 : NBBB=NBNN 926 : SEGINI WRK1,WRK2,WRK4 927 : I195=0 928 : I259=0 929 : C 930 : IG1=0 931 : C 932 : DO IB=1,NBELEM 933 : c 934 : c on cherche les coordonnees des noeuds de l element ib 935 : c 935 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 935 : CALL ZERO (REL,LRE,LRE) 936 : c 937 : c calcul des coordonnees locales de l'element 938 : c 938 : CALL JO2LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 939 : c 940 : c boucle sur les points de gauss 941 : c 942 : ISDJC=0 943 : DO IGAU=1,NBPGAU 944 : MPTVAL=IVAMAT 945 : DO IM=1,NMATT 946 : MELVAL=IVAL(IM) 947 : IGMN=MIN(IGAU,VELCHE(/1)) 948 : IBMN=MIN(IB,VELCHE(/2)) 949 : VALMAT(IM)=VELCHE(IGMN,IBMN) 950 : ENDDO 951 : C 952 : EPAIST=VALMAT(2) 953 : IF(EPAIST.EQ.0.D0)THEN 953 : CALL BJO2GN(IGAU,MFR,IFOUR,NIFOUR,XE,XEL,BPSS,SHPTOT, 954 : . SHPWRK,EPAIST,BGENE,DJAC,XDPGE,YDPGE,IERT) 955 : IF(IERT.NE.0) IG1=IB 956 : ENDIF 957 : C 957 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 958 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 959 : * 960 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 961 : IF(DJAC.EQ.0.) I259=IB 962 : DJAC=ABS(DJAC)*POIGAU(IGAU) 963 : * 964 : c valmat(1)=rho, valmat(2)=epai 965 : c /4 correspnnd en fait a diviser les matrices B par 2 966 : CCCCCCCCCCCC DJAC=DJAC*VALMAT(1)*VALMAT(2)/4.0D0 967 : DJAC=DJAC*VALMAT(1)*EPAIST/4.0D0 968 : c 969 : c cas axisymetrique : multiplication par le rayon de courbure 970 : c 971 : IF (IFOUR.EQ.0) THEN 972 : RAYON = 0.D0 973 : NUMSUP=NBNO/2 974 : DO IRAY=1,NUMSUP 975 : RAYON=RAYON + ( SHPTOT(1,IRAY,IGAU)*XE(1,IRAY) ) 976 : ENDDO 977 : DJAC=DJAC*RAYON 978 : ENDIF 979 : c 979 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 980 : ENDDO 981 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 982 : * SEGINI XMATRI 983 : * IMATTT(IB)=XMATRI 984 : c 985 : c remplissage de xmatri 986 : c 986 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 987 : * SEGDES XMATRI 988 : ENDDO 989 : C 990 : IF(IG1.NE.0) INTERR(1)=IG1 990 : IF(IG1.NE.0) CALL ERREUR (611) 991 : C 992 : IF(I195.NE.0) INTERR(1)=I195 992 : IF(I195.NE.0) CALL ERREUR(195) 993 : IF(I259.NE.0) INTERR(1)=I259 993 : IF(I259.NE.0) CALL ERREUR(259) 994 : SEGDES xMATRI 995 : SEGSUP WRK1,WRK2,WRK4,MVELCH 996 : GOTO 510 997 : C 998 : c_______________________________________________________________________ 999 : c 1000 : c secteur de calcul pour les elements joints jct3 en 2D cisaillement 1001 : c_______________________________________________________________________ 1002 : c 1003 : 168 CONTINUE 1004 : IF (ILUMP .EQ. 1) GOTO 99 1005 : NBNO=NBNN 1006 : NBBB=NBNN 1007 : SEGINI WRK1,WRK2,WRK4 1008 : I195=0 1009 : I259=0 1010 : DO IB=1,NBELEM 1011 : c 1012 : c on cherche les coordonnees des noeuds de l element ib 1013 : c 1013 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1013 : CALL ZERO (REL,LRE,LRE) 1014 : c 1015 : c calcul des coordonnees locales de l'element 1016 : c 1016 : CALL JT3LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 1017 : c 1018 : c boucle sur les points de gauss 1019 : c 1020 : ISDJC=0 1021 : DO IGAU=1,NBPGAU 1021 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1022 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 1023 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 1024 : IF(DJAC.EQ.0.) I259=IB 1025 : DJAC=ABS(DJAC)*POIGAU(IGAU) 1026 : MPTVAL=IVAMAT 1027 : IF (IVAL(1).NE.0) THEN 1028 : MELVAL=IVAL(1) 1029 : IGMN=MIN(IGAU,VELCHE(/1)) 1030 : IBMN=MIN(IB,VELCHE(/2)) 1031 : VALMAT(1)=VELCHE(IGMN,IBMN) 1032 : ELSE 1033 : VALMAT(1)=0.D0 1034 : ENDIF 1035 : DJAC=DJAC*VALMAT(1) 1036 : CCCCCCCCCC DJAC=DJAC/3.0D0 1037 : C Il faut diviser par 4, ce qui correspond plus exactement a diviser 1038 : C le B par 2... 1039 : DJAC=DJAC/4.0D0 1039 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 1040 : ENDDO 1041 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1042 : * SEGINI XMATRI 1043 : * IMATTT(IB)=XMATRI 1044 : c 1045 : c remplissage de xmatri 1046 : c 1046 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1047 : SEGDES XMATRI 1048 : ENDDO 1049 : IF(I195.NE.0) INTERR(1)=I195 1049 : IF(I195.NE.0) CALL ERREUR(195) 1050 : IF(I259.NE.0) INTERR(1)=I259 1050 : IF(I259.NE.0) CALL ERREUR(259) 1051 : SEGDES xMATRI 1052 : SEGSUP WRK1,WRK2,WRK4,MVELCH 1053 : GOTO 510 1054 : c_______________________________________________________________________ 1055 : c 1056 : c secteur de calcul pour les elements joints jgt3 generalise 1057 : c_______________________________________________________________________ 1058 : c 1059 : 171 CONTINUE 1060 : IF (ILUMP .EQ. 1) GOTO 99 1061 : NBNO=NBNN 1062 : NBBB=NBNN 1063 : SEGINI WRK1,WRK2,WRK4 1064 : I195=0 1065 : I259=0 1066 : C 1067 : IG1=0 1068 : C 1069 : DO IB=1,NBELEM 1070 : c 1071 : c on cherche les coordonnees des noeuds de l element ib 1072 : c 1072 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1072 : CALL ZERO (REL,LRE,LRE) 1073 : c 1074 : c calcul des coordonnees locales de l'element 1075 : c 1075 : CALL JT3LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 1076 : c 1077 : c boucle sur les points de gauss 1078 : c 1079 : ISDJC=0 1080 : DO IGAU=1,NBPGAU 1081 : MPTVAL=IVAMAT 1082 : DO IM=1,NMATT 1083 : MELVAL=IVAL(IM) 1084 : IGMN=MIN(IGAU,VELCHE(/1)) 1085 : IBMN=MIN(IB,VELCHE(/2)) 1086 : VALMAT(IM)=VELCHE(IGMN,IBMN) 1087 : ENDDO 1088 : C 1089 : EPAIST=VALMAT(2) 1090 : IF(EPAIST.EQ.0.D0)THEN 1090 : CALL BJT3G(IGAU,MFR,IFOUR,NIFOUR,XE,XEL,BPSS,SHPTOT, 1091 : . SHPWRK,EPAIST,BGENE,DJAC,IERT) 1092 : IF(IERT.NE.0) IG1=IB 1093 : ENDIF 1094 : C 1094 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1095 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 1096 : * 1097 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 1098 : IF(DJAC.EQ.0.) I259=IB 1099 : DJAC=ABS(DJAC)*POIGAU(IGAU) 1100 : * 1101 : c valmat(1)=rho, valmat(2)=epai 1102 : c /4 correspnnd en fait a diviser les matrices B par 2 1103 : CCCCCCCCCCCC DJAC=DJAC*VALMAT(1)*VALMAT(2)/4.0D0 1104 : DJAC=DJAC*VALMAT(1)*EPAIST/4.0D0 1105 : * 1105 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 1106 : ENDDO 1107 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1108 : * SEGINI XMATRI 1109 : * IMATTT(IB)=XMATRI 1110 : c 1111 : c remplissage de xmatri 1112 : c 1112 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1113 : * SEGDES XMATRI 1114 : ENDDO 1115 : C 1116 : IF(IG1.NE.0) INTERR(1)=IG1 1116 : IF(IG1.NE.0) CALL ERREUR (611) 1117 : C 1118 : IF(I195.NE.0) INTERR(1)=I195 1118 : IF(I195.NE.0) CALL ERREUR(195) 1119 : IF(I259.NE.0) INTERR(1)=I259 1119 : IF(I259.NE.0) CALL ERREUR(259) 1120 : SEGDES xMATRI 1121 : SEGSUP WRK1,WRK2,WRK4,MVELCH 1122 : GOTO 510 1123 : C 1124 : c_______________________________________________________________________ 1125 : c 1126 : c secteur de calcul pour les elements joints jci4 en 2D cisaillement 1127 : c_______________________________________________________________________ 1128 : c 1129 : 169 CONTINUE 1130 : IF (ILUMP .EQ. 1) GOTO 99 1131 : NBNO=NBNN 1132 : NBBB=NBNN 1133 : SEGINI WRK1,WRK2,WRK4 1134 : I195=0 1135 : I259=0 1136 : DO IB=1,NBELEM 1137 : c 1138 : c on cherche les coordonnees des noeuds de l element ib 1139 : c 1139 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1139 : CALL ZERO (REL,LRE,LRE) 1140 : c 1141 : c calcul des coordonnees locales de l'element 1142 : c 1142 : CALL JO4LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 1143 : c 1144 : c boucle sur les points de gauss 1145 : c 1146 : ISDJC=0 1147 : DO IGAU=1,NBPGAU 1147 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1148 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 1149 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 1150 : IF(DJAC.EQ.0.) I259=IB 1151 : DJAC=ABS(DJAC)*POIGAU(IGAU) 1152 : MPTVAL=IVAMAT 1153 : IF (IVAL(1).NE.0) THEN 1154 : MELVAL=IVAL(1) 1155 : IGMN=MIN(IGAU,VELCHE(/1)) 1156 : IBMN=MIN(IB,VELCHE(/2)) 1157 : VALMAT(1)=VELCHE(IGMN,IBMN) 1158 : ELSE 1159 : VALMAT(1)=0.D0 1160 : ENDIF 1161 : DJAC=DJAC*VALMAT(1) 1162 : CCCCCCCCCCC DJAC=DJAC/3.0D0 1163 : C Il faut diviser par 4, ce qui correspond plus exactement a diviser 1164 : C le B par 2... 1165 : DJAC=DJAC/4.0D0 1165 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 1166 : ENDDO 1167 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1168 : * SEGINI XMATRI 1169 : * IMATTT(IB)=XMATRI 1170 : c 1171 : c remplissage de xmatri 1172 : c 1172 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1173 : * SEGDES XMATRI 1174 : ENDDO 1175 : IF(I195.NE.0) INTERR(1)=I195 1175 : IF(I195.NE.0) CALL ERREUR(195) 1176 : IF(I259.NE.0) INTERR(1)=I259 1176 : IF(I259.NE.0) CALL ERREUR(259) 1177 : SEGDES xMATRI 1178 : SEGSUP WRK1,WRK2,WRK4,MVELCH 1179 : GOTO 510 1180 : c_______________________________________________________________________ 1181 : c 1182 : c secteur de calcul pour les elements joints jgi4 generalise 1183 : c_______________________________________________________________________ 1184 : c 1185 : 172 CONTINUE 1186 : IF (ILUMP .EQ. 1) GOTO 99 1187 : NBNO=NBNN 1188 : NBBB=NBNN 1189 : SEGINI WRK1,WRK2,WRK4 1190 : I195=0 1191 : I259=0 1192 : C 1193 : IG1=0 1194 : C 1195 : DO IB=1,NBELEM 1196 : c 1197 : c on cherche les coordonnees des noeuds de l element ib 1198 : c 1198 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1198 : CALL ZERO (REL,LRE,LRE) 1199 : c 1200 : c calcul des coordonnees locales de l'element 1201 : c 1201 : CALL JO4LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) 1202 : c 1203 : c boucle sur les points de gauss 1204 : c 1205 : ISDJC=0 1206 : DO IGAU=1,NBPGAU 1207 : MPTVAL=IVAMAT 1208 : DO IM=1,NMATT 1209 : MELVAL=IVAL(IM) 1210 : IGMN=MIN(IGAU,VELCHE(/1)) 1211 : IBMN=MIN(IB,VELCHE(/2)) 1212 : VALMAT(IM)=VELCHE(IGMN,IBMN) 1213 : ENDDO 1214 : C 1215 : EPAIST=VALMAT(2) 1216 : IF(EPAIST.EQ.0.D0)THEN 1216 : CALL BJO4G(IGAU,XE,XEL,BPSS,SHPTOT,SHPWRK,EPAIST, 1217 : . BGENE,DJAC,IERT) 1218 : IF(IERT.NE.0) IG1=IB 1219 : ENDIF 1220 : C 1220 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1221 : 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 1222 : * 1223 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 1224 : IF(DJAC.EQ.0.) I259=IB 1225 : DJAC=ABS(DJAC)*POIGAU(IGAU) 1226 : * 1227 : c valmat(1)=rho, valmat(2)=epai 1228 : c /4 correspnnd en fait a diviser les matrices B par 2 1229 : CCCCCCCCCCCC DJAC=DJAC*VALMAT(1)*VALMAT(2)/4.0D0 1230 : DJAC=DJAC*VALMAT(1)*EPAIST/4.0D0 1231 : * 1231 : CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 1232 : ENDDO 1233 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1234 : * SEGINI XMATRI 1235 : * IMATTT(IB)=XMATRI 1236 : c 1237 : c remplissage de xmatri 1238 : c 1238 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1239 : * SEGDES XMATRI 1240 : ENDDO 1241 : C 1242 : IF(IG1.NE.0) INTERR(1)=IG1 1242 : IF(IG1.NE.0) CALL ERREUR (611) 1243 : C 1244 : IF(I195.NE.0) INTERR(1)=I195 1244 : IF(I195.NE.0) CALL ERREUR(195) 1245 : IF(I259.NE.0) INTERR(1)=I259 1245 : IF(I259.NE.0) CALL ERREUR(259) 1246 : SEGDES xMATRI 1247 : SEGSUP WRK1,WRK2,WRK4,MVELCH 1248 : GOTO 510 1249 : C 1250 : c_______________________________________________________________________ 1251 : c 1252 : c secteur de calcul pour les elements homogeneises 1253 : c (liquide solide) trih 1254 : c_______________________________________________________________________ 1255 : c 1256 : 92 CONTINUE 1257 : IF (ILUMP .EQ. 1) GOTO 99 1258 : NBNO=NBNN 1259 : NBBB=NBNN 1260 : SEGINI WRK1,WRK2 1261 : I195=0 1262 : DO 3092 IB=1,NBELEM 1263 : c 1264 : c on cherche les coordonnees des noeuds de l element ib 1265 : c 1265 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1265 : CALL ZERO (REL,LRE,LRE) 1266 : c 1267 : c on cherche les caracteristiques du materiau de l element ib 1268 : c 1269 : MPTVAL=IVAMAT 1270 : DO 8092 IM=1,NMATT 1271 : MELVAL=IVAL(IM) 1272 : IBMN=MIN(IB,VELCHE(/2)) 1273 : VALMAT(IM)=VELCHE(1,IBMN) 1274 : 8092 CONTINUE 1275 : B11 =VALMAT(1) 1276 : B22 =VALMAT(2) 1277 : B12 =VALMAT(3) 1278 : RHOF =VALMAT(4) 1279 : RHOS =VALMAT(5) 1280 : C =VALMAT(6) 1281 : RHOREF=VALMAT(7) 1282 : CREF =VALMAT(8) 1283 : RLCAR =VALMAT(9) 1284 : c 1285 : c on cherche les caracteristiques geometriques de l element ib 1286 : c 1287 : MPTVAL=IVACAR 1288 : IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN 1289 : MELVAL=IVAL(1) 1290 : IBMN=MIN(IB,VELCHE(/2)) 1291 : SECT=VELCHE(1,IBMN) 1292 : MELVAL=IVAL(2) 1293 : IBMN=MIN(IB,VELCHE(/2)) 1294 : SCEL=VELCHE(1,IBMN) 1295 : MELVAL=IVAL(3) 1296 : IBMN=MIN(IB,VELCHE(/2)) 1297 : SFLU=VELCHE(1,IBMN) 1298 : MELVAL=IVAL(4) 1299 : IBMN=MIN(IB,VELCHE(/2)) 1300 : EPS =VELCHE(1,IBMN) 1301 : ELSE 1302 : SECT=1.D0 1303 : MELVAL=IVAL(1) 1304 : IBMN=MIN(IB,VELCHE(/2)) 1305 : SCEL=VELCHE(1,IBMN) 1306 : MELVAL=IVAL(2) 1307 : IBMN=MIN(IB,VELCHE(/2)) 1308 : SFLU=VELCHE(1,IBMN) 1309 : MELVAL=IVAL(3) 1310 : IBMN=MIN(IB,VELCHE(/2)) 1311 : EPS =VELCHE(1,IBMN) 1312 : MELVAL=IVAL(4) 1313 : IBMN=MIN(IB,VELCHE(/2)) 1314 : F11 =VELCHE(1,IBMN) 1315 : MELVAL=IVAL(5) 1316 : IBMN=MIN(IB,VELCHE(/2)) 1317 : F12 =VELCHE(1,IBMN) 1318 : ENDIF 1319 : c 1320 : c calcul de la masse m0/eps**2 1321 : c 1322 : RHOSS=RHOS*SECT/(EPS*EPS) 1323 : c 1324 : c calcul des coefficients de normalisation 1325 : c 1326 : COEFPR=(RHOREF*CREF*CREF)/RLCAR 1327 : COEFPI=RHOREF*RLCAR 1328 : VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) 1329 : VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) 1330 : VKL23 = COEFPI/SCEL 1331 : VKL33 = 1.D0/SCEL 1332 : IF(IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN 1333 : VKL23 =COEFPI*0.5D0*(2.D0*SCEL-B11-B22)/SCEL 1334 : VKL33 =(RHOSS+RHOF*(SFLU-(B11+B22)/2.D0))/SCEL 1335 : c 1336 : c calcul des termes en pi*pi 1337 : c integration par nbpgau points de gauss 1338 : c 1339 : ISDJC=0 1340 : DO 4092 IGAU=1,NBPGAU 1341 : POIGA2=MINTE.POIGAU(IGAU) 1341 : CALL TRIHM1(IGAU,MELE,MFR,NBNO,XE,SHPTOT,SHPWRK, 1342 : # IFOUR,NHRM,B11,B22,SFLU,POIGA2,VKL22,LRE,REL,IRRT) 1343 : IF(IRRT.NE.1) GOTO 7092 1344 : c 1345 : c calcul du reste des termes de la matrice masse 1346 : c integration par nbpgau points de gauss 1347 : c 1347 : CALL TRIHM2(IGAU,MELE,MFR,NBNO,XE,MINTE.SHPTOT,SHPWRK 1348 : # ,IFOUR,NHRM,VKL12,VKL23,VKL33,POIGA2,ISDJC,LRE,REL,IRRT) 1349 : IF(IRRT.NE.1) GOTO 7092 1350 : 4092 CONTINUE 1351 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1352 : ELSE 1353 : c 1354 : c boucle sur les points de gauss 1355 : c cas plan 1356 : c 1357 : ISDJC=0 1358 : DO 6092 IGAU1=1,NBPGAU 1359 : POIGA1=MINTE.POIGAU(IGAU1) 1359 : CALL TRIHM3(IGAU1,MELE,NBNO,XE,SHPTOT,SHPWRK 1360 : # ,RHOSS,RHOF, 1361 : # B11,B22,B12,F11,F12,SFLU,SCEL,POIGA1,VKL12,VKL22, 1362 : # VKL23,VKL33,LRE,REL,ISDJC,IRRT) 1363 : IF(IRRT.NE.1) GOTO 7092 1364 : 6092 CONTINUE 1365 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1366 : ENDIF 1367 : * SEGINI XMATRI 1368 : * IMATTT(IB)=XMATRI 1369 : c 1370 : c remplissage de xmatri 1371 : c 1371 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1372 : * SEGDES XMATRI 1373 : 3092 CONTINUE 1374 : c 1375 : c impression d un eventuel message d erreur 1376 : c 1377 : IF(I195.NE.0) INTERR(1)=I195 1377 : IF(I195.NE.0) CALL ERREUR(195) 1378 : 7092 CONTINUE 1379 : IF(IRRT.EQ.0) THEN 1380 : MOTERR(1:4)=NOMTP(MELE) 1380 : CALL ERREUR(420) 1381 : ELSE 1382 : IF(IRRT.EQ.2) THEN 1383 : INTERR(1) = IB 1383 : CALL ERREUR(405) 1384 : ENDIF 1385 : ENDIF 1385 : IF(IRRT.EQ.3) CALL ERREUR(421) 1385 : IF(IRRT.EQ.4) CALL ERREUR(422) 1386 : SEGDES xMATRI 1387 : SEGSUP WRK1,WRK2,MVELCH 1388 : GOTO 510 1389 : c_______________________________________________________________________ 1390 : c 1391 : c secteur de calcul pour les elements de raccord 1392 : c liquide coque 4 noeuds - cas tridimensionnel 1393 : c_______________________________________________________________________ 1394 : c 1395 : 94 CONTINUE 1396 : IF (ILUMP .EQ. 1) GOTO 99 1397 : NBBB=NBNN 1398 : LW=IDIM 1399 : SEGINI WRK1,WRK3 1400 : DO 3094 IB=1,NBELEM 1401 : c 1402 : c on cherche les coordonnees de l element ib 1403 : c 1403 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1403 : CALL ZERO(REL,LRE,LRE) 1404 : c 1405 : c calcul des coefficients de normalisation 1406 : c 1407 : MPTVAL=IVAMAT 1408 : DO 5094 IM=1,NMATT 1409 : MELVAL=IVAL(IM) 1410 : IBMN=MIN(IB,VELCHE(/2)) 1411 : VALMAT(IM)=VELCHE(1,IBMN) 1412 : 5094 CONTINUE 1413 : RHOREF=VALMAT(1) 1414 : RLCAR = VALMAT(2) 1415 : c 1416 : CFPI= RHOREF*RLCAR 1417 : c 1418 : MPTVAL=IVACAR 1419 : DO 4094 IC=1,NCARR 1420 : IF (IVAL(IC).NE.0) THEN 1421 : MELVAL=IVAL(IC) 1422 : IBMN=MIN(IB,VELCHE(/2)) 1423 : WORK(IC)=VELCHE(1,IBMN) 1424 : ELSE 1425 : WORK(IC)=0.D0 1426 : ENDIF 1427 : 4094 CONTINUE 1428 : c 1428 : CALL LIC4MA(NBPGAU,IDIM,NBNN,NDDL,XE,CFPI,WORK,POIGAU, 1429 : 1 QSIGAU,ETAGAU,SHPTOT,REL,LRE,IER246) 1430 : IF(IER246.NE.0) THEN 1430 : CALL ERREUR(IER246) 1431 : SEGSUP xMATRI 1432 : SEGSUP WRK1,WRK3,MVELCH 1433 : GOTO 510 1434 : ENDIF 1435 : * SEGINI XMATRI 1436 : * IMATTT(IB)=XMATRI 1437 : c 1438 : c remplissage de xmatri 1439 : c 1439 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1440 : * SEGDES XMATRI 1441 : 3094 CONTINUE 1442 : SEGDES xMATRI 1443 : SEGSUP WRK1,WRK3,MVELCH 1444 : GOTO 510 1445 : c_______________________________________________________________________ 1446 : c 1447 : c secteur de calcul pour les elements homogeneises 1448 : c (liquide solide) quah 1449 : c_______________________________________________________________________ 1450 : c 1451 : 126 CONTINUE 1452 : c 1453 : IF (ILUMP .EQ. 1) GOTO 99 1454 : NBNO=NBNN 1455 : NBBB=NBNN 1456 : SEGINI WRK1,WRK2 1457 : I195=0 1458 : DO 3126 IB=1,NBELEM 1459 : c 1460 : c on cherche les coordonnees des noeuds de l element ib 1461 : c 1461 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1461 : CALL ZERO (REL,LRE,LRE) 1462 : c 1463 : c on cherche les caracteristiques du materiau de l element ib 1464 : c 1465 : MPTVAL=IVAMAT 1466 : DO 8126 IM=1,NMATT 1467 : MELVAL=IVAL(IM) 1468 : IBMN=MIN(IB,VELCHE(/2)) 1469 : VALMAT(IM)=VELCHE(1,IBMN) 1470 : 8126 CONTINUE 1471 : B11 =VALMAT(1) 1472 : B22 =VALMAT(2) 1473 : B12 =VALMAT(3) 1474 : RHOF =VALMAT(4) 1475 : RHOS =VALMAT(5) 1476 : C =VALMAT(6) 1477 : RHOREF=VALMAT(7) 1478 : CREF =VALMAT(8) 1479 : RLCAR =VALMAT(9) 1480 : c 1481 : c on cherche les caracteristiques geometriques de l element ib 1482 : c 1483 : MPTVAL=IVACAR 1484 : MELVAL=IVAL(4) 1485 : IBMN=MIN(IB,VELCHE(/2)) 1486 : SECT=VELCHE(1,IBMN) 1487 : c 1488 : MELVAL=IVAL(1) 1489 : IBMN=MIN(IB,VELCHE(/2)) 1490 : SCEL=VELCHE(1,IBMN) 1491 : c 1492 : MELVAL=IVAL(2) 1493 : IBMN=MIN(IB,VELCHE(/2)) 1494 : SFLU=VELCHE(1,IBMN) 1495 : c 1496 : MELVAL=IVAL(3) 1497 : IBMN=MIN(IB,VELCHE(/2)) 1498 : EPS =VELCHE(1,IBMN) 1499 : c 1500 : c 1501 : c calcul de la masse m0/eps**2 1502 : c 1503 : RHOSS=RHOS*SECT/(EPS*EPS) 1504 : c 1505 : c calcul des coefficients de normalisation 1506 : c 1507 : COEFPR=(RHOREF*CREF*CREF)/RLCAR 1508 : COEFPI=RHOREF*RLCAR 1509 : VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) 1510 : VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) 1511 : VKL23 =COEFPI*0.5D0*(2.D0*SCEL-B11-B22)/SCEL 1512 : VKL33 =(RHOSS+RHOF*(SFLU-(B11+B22)/2.D0))/SCEL 1513 : c 1514 : c 1515 : c calcul des termes en pi*pi 1516 : c integration par nbpgau points de gauss 1517 : c 1518 : ISDJC=0 1519 : DO 4126 IGAU=1,NBPGAU 1520 : POIGA2=MINTE.POIGAU(IGAU) 1520 : CALL QUAHM1(IGAU,MELE,MFR,NBNO,XE,SHPTOT,SHPWRK,IFOUR 1521 : # ,NHRM,B11,B22,SFLU,POIGA2,VKL22,LRE,REL,IRRT) 1522 : IF(IRRT.NE.1) GOTO 7126 1523 : c 1524 : c calcul du reste des termes de la matrice masse 1525 : c integration par nbpgau points de gauss 1526 : c 1526 : CALL QUAHM2(IGAU,MELE,MFR,NBNO,XE,MINTE.SHPTOT,SHPWRK 1527 : # ,IFOUR,NHRM,VKL12,VKL23,VKL33,POIGA2,ISDJC,LRE,REL,IRRT) 1528 : IF(IRRT.NE.1) GOTO 7126 1529 : 4126 CONTINUE 1530 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1531 : c 1532 : * SEGINI XMATRI 1533 : * IMATTT(IB)=XMATRI 1534 : c 1535 : c remplissage de xmatri 1536 : c 1536 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1537 : * SEGDES XMATRI 1538 : 3126 CONTINUE 1539 : c 1540 : c impression d un eventuel message d erreur 1541 : c 1542 : IF(I195.NE.0) INTERR(1)=I195 1542 : IF(I195.NE.0) CALL ERREUR(195) 1543 : 7126 CONTINUE 1544 : IF(IRRT.EQ.0) THEN 1545 : MOTERR(1:4)=NOMTP(MELE) 1545 : CALL ERREUR(420) 1546 : ELSE 1547 : IF(IRRT.EQ.2) THEN 1548 : INTERR(1) = IB 1548 : CALL ERREUR(405) 1549 : ENDIF 1550 : ENDIF 1550 : IF(IRRT.EQ.3) CALL ERREUR(421) 1550 : IF(IRRT.EQ.4) CALL ERREUR(422) 1551 : SEGDES xMATRI 1552 : SEGSUP WRK1,WRK2,MVELCH 1553 : GOTO 510 1554 : c 1555 : c_______________________________________________________________________ 1556 : c 1557 : c secteur de calcul pour les elements homogeneises 1558 : c (liquide solide) cubh 1559 : c_______________________________________________________________________ 1560 : c 1561 : 127 CONTINUE 1562 : IF (ILUMP .EQ. 1) GOTO 99 1563 : NBNO=NBNN 1564 : NBBB=NBNN 1565 : LW=IDIM 1566 : SEGINI WRK1,WRK2 1567 : I195=0 1568 : DO 3127 IB=1,NBELEM 1569 : c 1570 : c on cherche les coordonnees des noeuds de l element ib 1571 : c 1571 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1571 : CALL ZERO (REL,LRE,LRE) 1572 : c 1573 : c on cherche les caracteristiques du materiau de l element ib 1574 : c 1575 : MPTVAL=IVAMAT 1576 : DO 8127 IM=1,NMATT 1577 : MELVAL=IVAL(IM) 1578 : IBMN=MIN(IB,VELCHE(/2)) 1579 : VALMAT(IM)=VELCHE(1,IBMN) 1580 : 8127 CONTINUE 1581 : B11 =VALMAT(1) 1582 : B22 =VALMAT(2) 1583 : B12 =VALMAT(3) 1584 : RHOF =VALMAT(4) 1585 : RHOS =VALMAT(5) 1586 : C =VALMAT(6) 1587 : RHOREF=VALMAT(7) 1588 : CREF =VALMAT(8) 1589 : RLCAR =VALMAT(9) 1590 : c 1591 : c on cherche les caracteristiques geometriques de l element ib 1592 : c 1593 : MPTVAL=IVACAR 1594 : c 1595 : MELVAL=IVAL(1) 1596 : IBMN=MIN(IB,VELCHE(/2)) 1597 : SCEL=VELCHE(1,IBMN) 1598 : c 1599 : MELVAL=IVAL(2) 1600 : IBMN=MIN(IB,VELCHE(/2)) 1601 : SFLU=VELCHE(1,IBMN) 1602 : c 1603 : MELVAL=IVAL(3) 1604 : IBMN=MIN(IB,VELCHE(/2)) 1605 : EPS =VELCHE(1,IBMN) 1606 : c 1607 : MELVAL=IVAL(4) 1608 : IBMN=MIN(IB,VELCHE(/2)) 1609 : SECT =VELCHE(1,IBMN) 1610 : c 1611 : c calcul de la masse m0/eps**2 1612 : c 1613 : RHOSS=RHOS*SECT/(EPS*EPS) 1614 : c 1615 : c calcul des coefficients de normalisation 1616 : c 1617 : COEFPR=(RHOREF*CREF*CREF)/RLCAR 1618 : COEFPI=RHOREF*RLCAR 1619 : c 1620 : c 1621 : VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) 1622 : VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) 1623 : c 1624 : VKL23 = COEFPI/SCEL 1625 : VKL33 = 1.D0/SCEL 1626 : c 1627 : ISDJC=0 1628 : DO 6127 IGAU1=1,NBPGAU 1629 : POIGA1=MINTE.POIGAU(IGAU1) 1629 : CALL CUBHM1(IGAU1,MELE,NBNO,XE,SHPTOT,SHPWRK, 1630 : # RHOSS,RHOF,B11,B22,B12,SFLU,SCEL,POIGA1,VKL12,VKL22,VKL23, 1631 : # VKL33,LRE,REL,ISDJC,IRRT) 1632 : IF(IRRT.NE.1) GOTO 7127 1633 : 6127 CONTINUE 1634 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1635 : * SEGINI XMATRI 1636 : * IMATTT(IB)=XMATRI 1637 : c 1638 : c remplissage de xmatri 1639 : c 1639 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1640 : * SEGDES XMATRI 1641 : 3127 CONTINUE 1642 : c 1643 : c impression d un eventuel message d erreur 1644 : c 1645 : IF(I195.NE.0) INTERR(1)=I195 1645 : IF(I195.NE.0) CALL ERREUR(195) 1646 : 7127 CONTINUE 1647 : IF(IRRT.EQ.0) THEN 1648 : MOTERR(1:4)=NOMTP(MELE) 1648 : CALL ERREUR(420) 1649 : ELSE 1650 : IF(IRRT.EQ.2) THEN 1651 : INTERR(1) = IB 1651 : CALL ERREUR(405) 1652 : ENDIF 1653 : ENDIF 1653 : IF(IRRT.EQ.3) CALL ERREUR(421) 1653 : IF(IRRT.EQ.4) CALL ERREUR(422) 1654 : SEGDES xMATRI 1655 : SEGSUP WRK1,WRK2,MVELCH 1656 : GOTO 510 1657 : c_______________________________________________________________________ 1658 : c 1659 : c secteur de calcul pour les elements homogeneises 1660 : c (liquide solide) trh6 1661 : c_______________________________________________________________________ 1662 : c 1663 : 157 CONTINUE 1664 : IF (ILUMP .EQ. 1) GOTO 99 1665 : NBNO=NBNN 1666 : NBBB=NBNN 1667 : SEGINI WRK1,WRK2 1668 : I195=0 1669 : DO 3157 IB=1,NBELEM 1670 : c 1671 : c on cherche les coordonnees des noeuds de l element ib 1672 : c 1672 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 1672 : CALL ZERO (REL,LRE,LRE) 1673 : c 1674 : c on cherche les caracteristiques du materiau de l element ib 1675 : c 1676 : MPTVAL=IVAMAT 1677 : DO 8157 IM=1,NMATT 1678 : MELVAL=IVAL(IM) 1679 : IBMN=MIN(IB,VELCHE(/2)) 1680 : VALMAT(IM)=VELCHE(1,IBMN) 1681 : 8157 CONTINUE 1682 : B11 =VALMAT(1) 1683 : B22 =VALMAT(2) 1684 : B12 =VALMAT(3) 1685 : RHOF =VALMAT(4) 1686 : RHOS =VALMAT(5) 1687 : C =VALMAT(6) 1688 : RHOREF=VALMAT(7) 1689 : CREF =VALMAT(8) 1690 : RLCAR =VALMAT(9) 1691 : E111 =VALMAT(10) 1692 : E112 =VALMAT(11) 1693 : E121 =VALMAT(12) 1694 : E122 =VALMAT(13) 1695 : E221 =VALMAT(14) 1696 : E222 =VALMAT(15) 1697 : c 1698 : c on cherche les caracteristiques geometriques de l element ib 1699 : c 1700 : MPTVAL=IVACAR 1701 : IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN 1702 : MELVAL=IVAL(1) 1703 : IBMN=MIN(IB,VELCHE(/2)) 1704 : SECT=VELCHE(1,IBMN) 1705 : MELVAL=IVAL(2) 1706 : IBMN=MIN(IB,VELCHE(/2)) 1707 : SCEL=VELCHE(1,IBMN) 1708 : MELVAL=IVAL(3) 1709 : IBMN=MIN(IB,VELCHE(/2)) 1710 : SFLU=VELCHE(1,IBMN) 1711 : MELVAL=IVAL(4) 1712 : IBMN=MIN(IB,VELCHE(/2)) 1713 : EPS =VELCHE(1,IBMN) 1714 : ELSE 1715 : SECT=1.D0 1716 : MELVAL=IVAL(1) 1717 : IBMN=MIN(IB,VELCHE(/2)) 1718 : SCEL=VELCHE(1,IBMN) 1719 : MELVAL=IVAL(2) 1720 : IBMN=MIN(IB,VELCHE(/2)) 1721 : SFLU=VELCHE(1,IBMN) 1722 : MELVAL=IVAL(3) 1723 : IBMN=MIN(IB,VELCHE(/2)) 1724 : EPS =VELCHE(1,IBMN) 1725 : MELVAL=IVAL(4) 1726 : IBMN=MIN(IB,VELCHE(/2)) 1727 : F11 =VELCHE(1,IBMN) 1728 : MELVAL=IVAL(5) 1729 : IBMN=MIN(IB,VELCHE(/2)) 1730 : F12 =VELCHE(1,IBMN) 1731 : ENDIF 1732 : c 1733 : c calcul de la masse m0/eps**2 1734 : c 1735 : RHOSS=RHOS*SECT/(EPS*EPS) 1736 : c 1737 : c calcul des coefficients de normalisation 1738 : c 1739 : COEFPR=(RHOREF*CREF*CREF)/RLCAR 1740 : COEFPI=RHOREF*RLCAR 1741 : VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) 1742 : VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) 1743 : VKL23 = COEFPI/SCEL 1744 : VKL33 = 1.D0/SCEL 1745 : VKL41 = EPS*EPS/2.D0/SCEL*(COEFPR*COEFPI) 1746 : VKL42 = EPS*EPS/2.D0/SCEL*COEFPI*COEFPI 1747 : VKL43 = EPS*EPS/2.D0/SCEL*COEFPI 1748 : VKL44 = EPS*EPS/2.D0/SCEL 1749 : c 1750 : c boucle sur les points de gauss 1751 : c cas plan 1752 : c 1753 : ISDJC=0 1754 : DO 6157 IGAU1=1,NBPGAU 1755 : POIGA1=MINTE.POIGAU(IGAU1) 1755 : CALL TRIHM31(IGAU1,MELE,NBNO,XE,SHPTOT,SHPWRK 1756 : # ,RHOSS,RHOF, 1757 : # B11,B22,B12,F11,F12,SFLU,SCEL,POIGA1,VKL12,VKL22, 1758 : # VKL23,VKL33,VKL42,VKL43,VKL44,E111,E112,E121,E122, 1759 : # E221,E222,LRE,REL,ISDJC,IRRT) 1760 : IF(IRRT.NE.1) GOTO 7092 1761 : 6157 CONTINUE 1762 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1763 : * SEGINI XMATRI 1764 : * IMATTT(IB)=XMATRI 1765 : c 1766 : c remplissage de xmatri 1767 : c 1767 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 1768 : * SEGDES XMATRI 1769 : 3157 CONTINUE 1770 : c 1771 : c impression d un eventuel message d erreur 1772 : c 1773 : IF(I195.NE.0) INTERR(1)=I195 1773 : IF(I195.NE.0) CALL ERREUR(195) 1774 : 7157 CONTINUE 1775 : IF(IRRT.EQ.0) THEN 1776 : MOTERR(1:4)=NOMTP(MELE) 1776 : CALL ERREUR(420) 1777 : ELSE 1778 : IF(IRRT.EQ.2) THEN 1779 : INTERR(1) = IB 1779 : CALL ERREUR(405) 1780 : ENDIF 1781 : ENDIF 1781 : IF(IRRT.EQ.3) CALL ERREUR(421) 1781 : IF(IRRT.EQ.4) CALL ERREUR(422) 1782 : SEGDES xMATRI 1783 : SEGSUP WRK1,WRK2,MVELCH 1784 : GOTO 510 1785 : c_______________________________________________________________________ 1786 : * 1787 : 99 CONTINUE 1788 : MOTERR(1:4)=NOMTP(MELE) 1789 : MOTERR(5:12)='MASSE4' 1789 : CALL ERREUR(86) 1790 : * 1791 : 510 CONTINUE 1792 : if (CMATE.EQ.'STATIQUE') then 1793 : mlmots = iinc 1794 : if (iinc.gt.0) segsup mlmots 1795 : mlmots = idua 1796 : if (idua.gt.0) segsup mlmots 1797 : endif 1798 : RETURN 1799 : END 1800 : 1801 : 1802 : 1803 : 1804 : 1805 : 1806 :
© Cast3M 2003 - All rights reserved.
Disclaimer