1 : C SIGMA1 SOURCE AM 14/02/06 21:15:05 7945 2 : 3 : SUBROUTINE SIGMA1(MATE,IMAT,IPMAIL,IPMINT,MELE,IELE, 4 : & IVADEP,NBPTEL,LRE,NSTRS,IVAMAT,NBGMAT,NELMAT,LHOOK,NMATT, 5 : & CMATE,MFR,NDEP,IPORE,IREPS2,NBPGAU,IVASTR,UZDPG,RYDPG,RXDPG, 6 : & IIPDPG) 7 : 8 : *---------------------------------------------------------------------* 9 : * _________________________ * 10 : * | | * 11 : * | calcul des contraintes | * 12 : * |_________________________| * 13 : * * 14 : * massif, poreux, joints poreux, incompressibles * 15 : * * 16 : *---------------------------------------------------------------------* 17 : * * 18 : * entrees : * 19 : * ________ * 20 : * * 21 : * mate numero du materiau * 22 : * imat (2 il y a une matrice de hooke,1 non ) * 23 : * ipmail pointeur sur un segment meleme * 24 : * ipmint pointeur sur un segment minte * 25 : * mele numero de l'element fini * 26 : * iele numero geometrique de l'element 27 : * nbpgau nombre de point d'integration pour la rigidite * 28 : * ivadep pointeur sur le chamelem de deplacements * 29 : * nbptel nombre de points par element * 30 : * lre nombre de ddl dans la matrice de rigidite * 31 : * nstrs nombre de composante de contraintes/deformations * 32 : * ivamat pointeur sur un segment mptval pour le materiau ou * 33 : * pour une matrice de hooke * 34 : * nbgmat taille maxi des melval du materiau (pt de gauss) * 35 : * nelmat taille maxi des melval du materiau (no d'element) * 36 : * lhook dimension de la matrice de hooke * 37 : * nmatt nombre de composante de materiau (imat=1) * 38 : * cmate nom du materiau * 39 : * mfr numero de la formulation de l'element fini * 40 : * ndep nombre de composantes de deplacements * 41 : * ipore nombre de fonctions de forme * 42 : * iresp2 flag pour indiquer si on veut les contraintes * 43 : * de piola-kirchhoff * 44 : * uzdpg = deformation au point nsdpge support de la * 45 : * rydpf = deformation plane generalisee * 46 : * rxdpg = * 47 : * * 48 : * sorties : * 49 : * ________ * 50 : * * 51 : * ivastr pointeur sur un segment mptval contenant les * 52 : * les melvals de contraints 53 : * * 54 : *---------------------------------------------------------------------* 55 : 56 : IMPLICIT INTEGER(I-N) 57 : IMPLICIT REAL*8(A-H,O-Z) 58 : 59 : -INC CCOPTIO 60 : -INC CCHAMP 61 : -INC CCREEL 62 : -INC CCGEOME 63 : -INC SMCHAML 64 : -INC SMINTE 65 : -INC SMELEME 66 : -INC SMCOORD 67 : -INC SMLREEL 68 : 69 : SEGMENT WRK1 70 : REAL*8 DDHOOK(LHOOK,LHOOK) ,XDDL(LRE) ,XSTRS(NSTRS) 71 : REAL*8 XE(3,NBBB) ,DDHOMU(LHOOK,LHOOK) 72 : ENDSEGMENT 73 : c 74 : SEGMENT WRK2 75 : REAL*8 SHPWRK(6,NBNO) ,BGENE(LHOOK,LRE) 76 : ENDSEGMENT 77 : c 78 : SEGMENT WRK3 79 : REAL*8 BPSS(3,3),XEL(3,NBBB) 80 : ENDSEGMENT 81 : c 82 : SEGMENT WRK5 83 : REAL*8 XGENE(NSTN,LRN),COBMA(LHOOK),XWRK(LHOOK) 84 : REAL*8 COBB(IDECAP),CPBB(IDECAP),KKBB(IDECAP,IDECAP) 85 : ENDSEGMENT 86 : c 87 : SEGMENT WRK8 88 : REAL*8 XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM),TXR(IDIM,IDIM) 89 : REAL*8 D1HO(LHOOK,LHOOK),ROTH(LHOOK,LHOOK) 90 : ENDSEGMENT 91 : c 92 : SEGMENT,MVELCH 93 : REAL*8 VALMAT(NV1) 94 : ENDSEGMENT 95 : c 96 : SEGMENT MPTVAL 97 : INTEGER IPOS(NS) ,NSOF(NS) 98 : INTEGER IVAL(NCOSOU) 99 : CHARACTER*16 TYVAL(NCOSOU) 100 : ENDSEGMENT 101 : * 102 : SEGMENT MTRACE 103 : REAL*8 TRACE(3,NBPTEL) 104 : ENDSEGMENT 105 : * 106 : CHARACTER*8 CMATE 107 : DIMENSION A(4,60),BB(3,60),UDPGE(3) 108 : LOGICAL BDPGE 109 : 110 : c Introduction du point autour duquel se fait le mouvement 111 : c de la section en defo plane generalisee 112 : c Pas de rotation en 1D 113 : BDPGE=.FALSE. 114 : NDPGE=0 115 : XDPGE=XZero 116 : YDPGE=XZero 117 : IF (IFOUR.EQ.-3) THEN 118 : BDPGE=.TRUE. 119 : NDPGE=3 120 : UDPGE(1)=UZDPG 121 : UDPGE(2)=RYDPG 122 : UDPGE(3)=RXDPG 123 : SEGACT,MCOORD 124 : IREF=(IIPDPG-1)*(IDIM+1) 125 : XDPGE=XCOOR(IREF+1) 126 : YDPGE=XCOOR(IREF+2) 127 : ELSE IF (IDIM.EQ.1) THEN 128 : IF ((IFOUR.GE.7 .AND. IFOUR.LE.10) .OR. IFOUR.EQ.14) THEN 129 : BDPGE=.TRUE. 130 : NDPGE=1 131 : UDPGE(1)=UZDPG 132 : ELSE IF (IFOUR.EQ.11) THEN 133 : BDPGE=.TRUE. 134 : NDPGE=2 135 : UDPGE(1)=UZDPG 136 : UDPGE(2)=RXDPG 137 : ENDIF 138 : ENDIF 139 : * 140 : MELEME=IPMAIL 141 : NBNN=NUM(/1) 142 : NBELEM=NUM(/2) 143 : * 144 : IDECAP=0 145 : NHRM=NIFOUR 146 : * 147 : NV1=NMATT 148 : SEGINI,MVELCH 149 : * 150 : MINTE=IPMINT 151 : * 152 : NBBB=NBNN 153 : SEGINI WRK1 154 : * 155 : IRTD=1 156 : c_______________________________________________________________________ 157 : c 158 : c numero des etiquettes : 159 : c etiquettes de 1 a 98 pour traitement specifique a l element 160 : c dans la zone specifique a chaque element commencant par : 161 : c 5 continue 162 : c element 5 etiquettes 1005 2005 3005 4005 ... 163 : c 44 continue 164 : c element 44 etiquettes 1044 2044 3044 4044 ... 165 : c_______________________________________________________________________ 166 : c 167 : GOTO (99,99,99, 4,99, 4,99, 4,99, 4,99,99,99, 4, 4, 4, 4,99,99,99, 168 : 1 99,99, 4, 4, 4, 4,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 169 : 2 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 170 : 3 99,99,99,99,99,99,99,99, 4, 4, 4, 4, 4,99,99, 4,99,99,79,79, 171 : 4 79,79,79,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 172 : 5 99,99,99,99,99,99,99,80,80,80, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 173 : 6 4, 4),MELE 174 : * 175 : IF (MELE.EQ.183.OR.MELE.EQ.184.OR. 176 : . MELE.EQ.193.OR.MELE.EQ.194) GOTO 4 177 : IF (MELE.GE.173.AND.MELE.LE.182) GOTO 173 178 : * 179 : GOTO 99 180 : 181 : c_______________________________________________________________________ 182 : c 183 : c elements massifs et elements incompressibles 184 : c_______________________________________________________________________ 185 : c 186 : 4 CONTINUE 187 : c 188 : c Cas non isotropes : 189 : c Recuperation des fonctions de forme et leurs derivees au centre de 190 : c l'element pour le calcul des axes locaux 191 : c 192 : IPMIN2 = 0 193 : IF ( (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) .AND. IMAT.EQ.1 ) THEN 194 : NLG=NUMGEO(MELE) 194 : CALL RESHPT(1,NBNN,NLG,MELE,0,IPMIN2,IRT1) 195 : MINTE2=IPMIN2 196 : SEGACT,MINTE2 197 : SEGINI,WRK8 198 : ENDIF 199 : c 200 : NBNO=NBNN 201 : SEGINI WRK2 202 : IF (IREPS2.EQ.1) SEGINI MTRACE 203 : c 204 : c* NDDD=NDEP 205 : c* IF (IFOUR.EQ.-3) NDDD=NDEP-3 206 : NDDD=NDEP-NDPGE 207 : c 208 : DO 3004 IB=1,NBELEM 209 : c 210 : c on cherche les deplacements 211 : c 212 : MPTVAL=IVADEP 213 : IE=1 214 : DO 4004 IGAU=1,NBNN 215 : DO 4004 ICOMP=1,NDDD 216 : MELVAL=IVAL(ICOMP) 217 : IGMN=MIN(IGAU,VELCHE(/1)) 218 : IBMN=MIN(IB ,VELCHE(/2)) 219 : XDDL(IE)=VELCHE(IGMN,IBMN) 220 : IE=IE+1 221 : 4004 CONTINUE 222 : IF (BDPGE) THEN 223 : DO i=1,NDPGE 224 : XDDL(IE)=UDPGE(i) 225 : IE=IE+1 226 : ENDDO 227 : ENDIF 228 : c 229 : c on cherche les coordonnees des noeuds de l element ib 230 : c 230 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 231 : c 232 : c calcul des axes locaux dans le cas des materiaux orthotropes, 233 : c anisotropes et unidirectionnel 234 : c 235 : C* IF ((MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) .AND. IMAT.EQ.1) THEN 236 : IF (IPMIN2.NE.0) THEN 237 : NBSH=MINTE2.SHPTOT(/2) 237 : CALL RLOCAL(XE,MINTE2.SHPTOT,NBSH,NBNN,TXR) 238 : IF (nbsh.EQ.-1) THEN 238 : CALL ERREUR(525) 239 : GOTO 9904 240 : ENDIF 241 : ENDIF 242 : c 243 : c calcul des coeff de modification de la matrice b-barre (incompres) 244 : c 245 : IF (MFR.EQ.31) THEN 245 : CALL BBCALC(XE,MELE,NBNN,IDIM,NBPGAU, 246 : 1 POIGAU,QSIGAU,ETAGAU,DZEGAU,NSTRS, 247 : 2 LRE,IFOUR,A,BB,NHRM,SHPTOT,SHPWRK,XDPGE,YDPGE) 248 : ENDIF 249 : c 250 : c boucle sur les points de gauss 251 : c 252 : ISDJC=0 253 : c 254 : DO 5004 IGAU=1,NBPTEL 255 : c 255 : CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 256 : 1 MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE, 257 : 2 SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 258 : c 259 : IF (DJAC.EQ.0.D0) THEN 260 : INTERR(1)=IB 260 : CALL ERREUR(259) 261 : GOTO 9904 262 : ELSE IF (DJAC.LT.0.D0) THEN 263 : ISDJC=ISDJC+1 264 : ENDIF 265 : 266 : C En cas d'elements incompressibles : BGENE selon la methode B-BARRE 267 : IF (MFR.EQ.31) THEN 267 : CALL BBAR(IGAU,NBPGAU, POIGAU,QSIGAU,ETAGAU,DZEGAU, 268 : & MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BB,BGENE) 269 : ENDIF 270 : c 271 : c on cherche les matrices de Hooke 272 : c 273 : MPTVAL=IVAMAT 274 : IF (IMAT.EQ.2) THEN 275 : MELVAL=IVAL(1) 276 : IBMN=MIN(IB ,IELCHE(/2)) 277 : IGMN=MIN(IGAU,IELCHE(/1)) 278 : MLREEL=IELCHE(IGMN,IBMN) 279 : IF (IGAU.LE.NBGMAT.AND.(IB.LE.NELMAT.OR.NBGMAT.GT.1)) THEN 280 : SEGACT MLREEL 280 : CALL DOHOOO(PROG,LHOOK,DDHOOK) 281 : SEGDES MLREEL 282 : ENDIF 283 : ELSE IF (IMAT.EQ.1) THEN 284 : DO 9004 IM=1,NMATT 285 : IF (IVAL(IM).NE.0) THEN 286 : MELVAL=IVAL(IM) 287 : IBMN=MIN(IB ,VELCHE(/2)) 288 : IGMN=MIN(IGAU,VELCHE(/1)) 289 : VALMAT(IM)=VELCHE(IGMN,IBMN) 290 : ELSE 291 : VALMAT(IM)=0.D0 292 : ENDIF 293 : 9004 CONTINUE 294 : c 295 : IF (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 296 : IF (IGAU.LE.NBGMAT) 296 : 1 CALL DOHMAO(VALMAT,CMATE,IFOUR,IDIM,TXR,XLOC,XGLOB,D1HO, 297 : 2 ROTH,DDHOOK,LHOOK,1,IRTD) 298 : ELSE 299 : IF (IGAU.LE.NBGMAT.AND.(IB.LE.NELMAT.OR.NBGMAT.GT.1)) 299 : 1 CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD) 300 : ENDIF 301 : ENDIF 302 : c 302 : CALL DBST(BGENE,DDHOOK,XDDL,LRE,NSTRS,XSTRS) 303 : c 304 : c calcul des eps 2 305 : c 306 : IF (IREPS2.EQ.1) 306 : 1 CALL DBST2(SHPWRK,DDHOOK,XDDL,XE,NBNO,IFOUR,NSTRS,XSTRS, 307 : 2 TRACE,IGAU,XDPGE,YDPGE,UDPGE,NHRM) 308 : c 309 : c remplissage du segment contenant les contraintes 310 : c 311 : MPTVAL=IVASTR 312 : DO 7004 ICOMP=1,NSTRS 313 : MELVAL=IVAL(ICOMP) 314 : IBMN=MIN(IB,VELCHE(/2)) 315 : VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 316 : 7004 CONTINUE 317 : c 318 : 5004 CONTINUE 319 : c 320 : IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN 321 : INTERR(1)=IB 321 : CALL ERREUR(195) 322 : GOTO 9904 323 : ENDIF 324 : c 325 : c Correction sur la partie quadratique de la contrainte dans le cas 326 : c des elements incompressibles 327 : c 328 : IF (IREPS2.EQ.1) THEN 329 : IF (MFR.EQ.31) THEN 329 : CALL DBBST2(TRACE,NBPTEL,IFOUR,MELE,POIGAU,QSIGAU, 330 : & ETAGAU,DZEGAU,SHPTOT,NBNO,SHPWRK,XE) 331 : L=2 332 : IF (IDIM.EQ.3 .OR. IFOUR .EQ. 0) L=3 333 : DO 5005 ICOMP=1,L 334 : MELVAL=IVAL(ICOMP) 335 : IBMN=MIN(IB ,VELCHE(/2)) 336 : DO 5006 IGAU=1,NBPTEL 337 : VELCHE(IGAU,IBMN)=VELCHE(IGAU,IBMN)+TRACE(1,IGAU) 338 : 5006 CONTINUE 339 : 5005 CONTINUE 340 : IF (L.EQ.2) THEN 341 : MELVAL=IVAL(3) 342 : IBMN=MIN(IB ,VELCHE(/2)) 343 : DO 5007 IGAU=1,NBPTEL 344 : VELCHE(IGAU,IBMN) = VELCHE(IGAU,IBMN) 345 : & + (TRACE(1,IGAU)/TRACE(2,IGAU)*TRACE(3,IGAU)) 346 : 5007 CONTINUE 347 : ENDIF 348 : ENDIF 349 : ENDIF 350 : 351 : 3004 CONTINUE 352 : c 353 : IF (IRTD.EQ.0) THEN 354 : MOTERR(1:8)=CMATE 355 : MOTERR(9:12)=NOMFR(MFR/2+1) 356 : INTERR(1)=IFOUR 356 : CALL ERREUR(81) 357 : ENDIF 358 : c 359 : 9904 CONTINUE 360 : SEGSUP WRK2 361 : IF (IREPS2.EQ.1) SEGSUP,MTRACE 362 : C* IF ((MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) .AND. IMAT.EQ.1) THEN 363 : IF (IPMIN2.NE.0) THEN 364 : SEGDES MINTE2 365 : SEGSUP WRK8 366 : ENDIF 367 : GOTO 510 368 : 369 : c____________________________________________________________________ 370 : c 371 : c milieux poreux 372 : c____________________________________________________________________ 373 : c 374 : 79 CONTINUE 375 : c 376 : c Ces cas ne sont pas prevus actuellement ! 377 : IF ( IMAT.EQ.2 .OR. 378 : & (IMAT.EQ.1.AND.(MATE.LT.1.OR.MATE.GT.4)) 379 : & ) GOTO 99 380 : 381 : c pour ces elements nbbb = nombre de noeuds 382 : c nbno = nombre de fonctions de forme 383 : c 384 : NBNO=IPORE 385 : NSTN=1 386 : LPP=0 387 : c***************** AM 08/01/01 388 : c* NSTMU=2 389 : c* IF (IFOUR.GE.0) NSTMU=3 390 : NSTMU=3 391 : LRN=NBNO-NBBB 392 : LRB=LRE-LRN 393 : c 394 : c Cas non isotropes : 395 : c recuperation des fonctions de forme et leurs derivees 396 : c au centre de l'element pour le calcul des axes locaux 397 : c 398 : IPMIN2 = 0 399 : IF ((MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) .AND. IMAT.EQ.1) THEN 399 : CALL RESHPT(1,NBNO,IELE,MELE,0,IPMIN2,IRT1) 400 : MINTE2=IPMIN2 401 : SEGACT MINTE2 402 : SEGINI WRK8 403 : NSTMU=LHOOK 404 : ENDIF 405 : c 406 : SEGINI WRK2,WRK5 407 : c Segment MTRACE initialise ici, necessaire mais inutilise 408 : IF (IREPS2.EQ.1) SEGINI MTRACE 409 : c 410 : I19 =0 411 : c 412 : DO 3079 IB=1,NBELEM 413 : c 414 : c on cherche d'abord les deplacements 415 : c 416 : MPTVAL=IVADEP 417 : IE=1 418 : DO 4079 IGAU=1,NBNN 419 : DO 4079 ICOMP=1,NDEP-1 420 : MELVAL=IVAL(ICOMP) 421 : IGMN=MIN(IGAU,VELCHE(/1)) 422 : IBMN=MIN(IB ,VELCHE(/2)) 423 : XDDL(IE)=VELCHE(IGMN,IBMN) 424 : IE=IE+1 425 : 4079 CONTINUE 426 : c 427 : c puis les pressions 428 : c 429 : MELVAL=IVAL(NDEP) 430 : DO 40791 IGAU=1,LRN 431 : IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1) 432 : IBMN=MIN(IB ,VELCHE(/2)) 433 : IGMN=MIN(IGAUSO,VELCHE(/1)) 434 : XDDL(IE)=VELCHE(IGMN,IBMN) 435 : IE=IE+1 436 : 40791 CONTINUE 437 : c 438 : c on cherche les coordonnees des noeuds de l element ib 439 : c 439 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 440 : c 441 : c calcul des axes locaux dans le cas des materiaux orthotropes, 442 : c anisotropes et unidirectionnels 443 : c 444 : C* IF ((MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) .AND. IMAT.EQ.1) THEN 445 : IF (IPMIN2.NE.0) THEN 446 : NBSH=MINTE2.SHPTOT(/2) 446 : CALL RLOCAL (XE,MINTE2.SHPTOT,NBSH,NBNN,TXR) 447 : if (nbsh.eq.-1) then 447 : call erreur(525) 448 : GOTO 9979 449 : endif 450 : ENDIF 451 : c 452 : c boucle sur les points de gauss 453 : c 454 : ISDJC=0 454 : CALL ZERO(COBMA,LHOOK,1) 455 : DO 5079 IGAU=1,NBPTEL 456 : c 456 : CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,NHRM, 457 : . 1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,1) 458 : c 459 : IF (DJAC.EQ.0.D0) THEN 460 : INTERR(1)=IB 460 : CALL ERREUR(259) 461 : GOTO 9979 462 : ELSE IF (DJAC.LT.0.D0) THEN 463 : ISDJC=ISDJC+1 464 : ENDIF 465 : c 466 : MPTVAL=IVAMAT 467 : C*D IF (IMAT.EQ.2) THEN 468 : C*D cas non prevu 469 : C*D GO TO 99 470 : C*D ELSE IF (IMAT.EQ.1) THEN 471 : DO 9079 IM=1,NMATT 472 : IF (IVAL(IM).NE.0) THEN 473 : MELVAL=IVAL(IM) 474 : IBMN=MIN(IB ,VELCHE(/2)) 475 : IGMN=MIN(IGAU,VELCHE(/1)) 476 : VALMAT(IM)=VELCHE(IGMN,IBMN) 477 : ELSE 478 : VALMAT(IM)=0.D0 479 : ENDIF 480 : 9079 CONTINUE 481 : IF (MATE.EQ.1) THEN 482 : IF (IGAU.LE.NBGMAT.AND.(IB.LE.NELMAT.OR.NBGMAT.GT.1)) 482 : . CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD) 483 : ELSE IF (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 484 : IF (IGAU.LE.NBGMAT) 484 : . CALL PORMAO(VALMAT,CMATE,IFOUR,IDIM,TXR,XLOC,XGLOB,D1HO, 485 : . ROTH,DDHOOK,LHOOK,COBMA,XMOB,1,IRTD) 486 : C*D ELSE 487 : C*D GOTO 99 488 : ENDIF 489 : C*D ENDIF 490 : c 490 : CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS) 491 : c 492 : c calcul des eps 2 493 : c 494 : IF (IREPS2.EQ.1) 494 : & CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,IGAU, 495 : & XDPGE,YDPGE,UDPGE,NHRM) 496 : * 497 : * contribution de epsi a msr0 498 : * 499 : IF (MATE.EQ.1) THEN 500 : C*D IF (IMAT.EQ.1) THEN 501 : DO 4879 I=1,NSTMU 502 : COBMA(I)=VALMAT(3) 503 : 4879 CONTINUE 504 : XMOB=VALMAT(4) 505 : C*D ELSE IF (IMAT.EQ.2) THEN 506 : C*D GO TO 99 507 : C*D ENDIF 508 : ENDIF 509 : * 510 : r_z=0. 511 : DO 4878 K=1,NSTMU 512 : r_z = r_z +COBMA(K)*XSTRS(K) 513 : 4878 CONTINUE 514 : XSTRS(NSTRS)=r_z 515 : DO 4877 KA=1,LHOOK 516 : XWRK(KA)=XSTRS(KA) 517 : 4877 CONTINUE 518 : * 519 : DO 4876 KA=1,LHOOK 520 : r_z =0.D0 521 : DO KB=1,LHOOK 522 : r_z = r_z + DDHOOK(KA,KB)*XWRK(KB) 523 : ENDDO 524 : XSTRS(KA)=r_z 525 : 4876 CONTINUE 526 : c 527 : c calcul de l'effet de la pression 528 : c 529 : IF (XMOB.EQ.0.D0) THEN 530 : UNSURM=0.D0 531 : ELSE 532 : UNSURM=1.D0 / XMOB 533 : ENDIF 534 : * 534 : CALL SIGPOR(COBMA,UNSURM,XGENE,NSTN,XDDL,IFOUR,NSTRS, 535 : . XSTRS,LRB,LRN,LPP,MELE,I19,COBB,KKBB,IDECAP) 536 : c 537 : c remplissage du segment contenant les contraintes 538 : c 539 : MPTVAL=IVASTR 540 : DO 7079 ICOMP=1,NSTRS 541 : MELVAL=IVAL(ICOMP) 542 : IBMN=MIN(IB ,VELCHE(/2)) 543 : VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 544 : 7079 CONTINUE 545 : c 546 : 5079 CONTINUE 547 : c 548 : IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN 549 : INTERR(1)=IB 549 : CALL ERREUR(195) 550 : GOTO 9979 551 : ENDIF 552 : c 553 : 3079 CONTINUE 554 : c 555 : IF (IRTD.EQ.0) THEN 556 : MOTERR(1:8)=CMATE 557 : MOTERR(9:12)=NOMFR(MFR/2+1) 558 : INTERR(1)=IFOUR 558 : CALL ERREUR(81) 559 : ENDIF 559 : IF (I19.NE.0) CALL ERREUR(19) 560 : c 561 : 9979 CONTINUE 562 : SEGSUP WRK2,WRK5 563 : IF (IREPS2.EQ.1) SEGSUP MTRACE 564 : C* IF ((MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) .AND. IMAT.EQ.1) THEN 565 : IF (IPMIN2.NE.0) THEN 566 : SEGDES MINTE2 567 : SEGSUP WRK8 568 : ENDIF 569 : GOTO 510 570 : 571 : c____________________________________________________________________ 572 : c 573 : c milieux poreux SUITE 574 : c____________________________________________________________________ 575 : c 576 : 173 CONTINUE 577 : c 578 : c Ces cas ne sont pas prevus actuellement ! 579 : IF ( IMAT.EQ.2 .OR. (IMAT.EQ.1.AND.MATE.NE.1) ) GOTO 99 580 : c 581 : c pour ces elements nbbb = nombre de noeuds 582 : c nbno = nombre de fonctions de forme 583 : c 584 : IF (MELE.GE.173.AND.MELE.LE.177) THEN 585 : IDECAP = 2 586 : ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN 587 : IDECAP = 3 588 : ENDIF 589 : * 590 : NBNO=IPORE 591 : NSTN=IDECAP 592 : 593 : LPP=NBNO-NBBB 594 : LRN=IDECAP*LPP 595 : LRB=LRE-LRN 596 : 597 : UNSURM = 0. 598 : 599 : SEGINI WRK2,WRK5 600 : c Segment MTRACE initialise ici (necessaire mais inutilise) 601 : IF (IREPS2.EQ.1) SEGINI MTRACE 602 : 603 : I19 =0 604 : c 605 : DO 3173 IB=1,NBELEM 606 : c 607 : c on cherche d'abord les deplacements 608 : c 609 : IE=1 610 : MPTVAL=IVADEP 611 : DO 4173 IGAU=1,NBNN 612 : DO 4173 ICOMP=1,NDEP-IDECAP 613 : MELVAL=IVAL(ICOMP) 614 : IGMN=MIN(IGAU,VELCHE(/1)) 615 : IBMN=MIN(IB ,VELCHE(/2)) 616 : XDDL(IE)=VELCHE(IGMN,IBMN) 617 : IE=IE+1 618 : 4173 CONTINUE 619 : c 620 : c puis les pressions 621 : c 622 : DO 4473 IPR = 1,IDECAP 623 : MELVAL=IVAL(NDEP-IDECAP+IPR) 624 : DO 4273 IGAU=1,LPP 625 : IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1) 626 : IBMN=MIN(IB ,VELCHE(/2)) 627 : IGMN=MIN(IGAUSO,VELCHE(/1)) 628 : XDDL(IE)=VELCHE(IGMN,IBMN) 629 : IE=IE+1 630 : 4273 CONTINUE 631 : 4473 CONTINUE 632 : c 633 : c on cherche les coordonnees des noeuds de l element ib 634 : c 634 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 635 : c 636 : c boucle sur les points de gauss 637 : c 638 : ISDJC=0 638 : CALL ZERO (COBB,IDECAP,1) 638 : CALL ZERO (CPBB,IDECAP,1) 638 : CALL ZERO (KKBB,IDECAP,IDECAP) 639 : c 640 : DO 5173 IGAU=1,NBPTEL 641 : c 641 : CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,NHRM, 642 : & 1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,1) 643 : c 644 : IF (DJAC.EQ.0.D0) THEN 645 : INTERR(1)=IB 645 : CALL ERREUR(259) 646 : GOTO 9973 647 : ELSE IF (DJAC.LT.0.D0) THEN 648 : ISDJC=ISDJC+1 649 : ENDIF 650 : c 651 : MPTVAL=IVAMAT 652 : C*D IF (IMAT.EQ.2) THEN 653 : C*D GO TO 99 654 : C*D ELSE IF (IMAT.EQ.1) THEN 655 : DO 6173 IM=1,NMATT 656 : IF (IVAL(IM).NE.0) THEN 657 : MELVAL=IVAL(IM) 658 : IBMN=MIN(IB ,VELCHE(/2)) 659 : IGMN=MIN(IGAU,VELCHE(/1)) 660 : VALMAT(IM)=VELCHE(IGMN,IBMN) 661 : ELSE 662 : VALMAT(IM)=0.D0 663 : ENDIF 664 : 6173 CONTINUE 665 : C*D IF (MATE.EQ.1) THEN 666 : IF (IGAU.LE.NBGMAT.AND.(IB.LE.NELMAT.OR.NBGMAT.GT.1)) 666 : . CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD) 667 : C*D ELSE 668 : C*D GOTO 99 669 : C*D ENDIF 670 : C*D ENDIF 671 : c 672 : IF (MFR.EQ.57) THEN 673 : COBB(1) = VALMAT(3) 674 : COBB(2) = VALMAT(4) 675 : CPBB(1) = VALMAT(5) 676 : CPBB(2) = VALMAT(6) 677 : KKBB(1,1)= VALMAT(7) 678 : KKBB(1,2)= VALMAT(8) 679 : KKBB(2,1)= VALMAT(9) 680 : KKBB(2,2)= VALMAT(10) 681 : ELSE IF(MFR.EQ.59) THEN 682 : COBB(1) = VALMAT(3) 683 : COBB(2) = VALMAT(4) 684 : COBB(3) = VALMAT(5) 685 : CPBB(1) = VALMAT(6) 686 : CPBB(2) = VALMAT(7) 687 : CPBB(3) = VALMAT(8) 688 : KKBB(1,1)= VALMAT(9) 689 : KKBB(1,2)= VALMAT(10) 690 : KKBB(1,3)= VALMAT(11) 691 : KKBB(2,1)= VALMAT(12) 692 : KKBB(2,2)= VALMAT(13) 693 : KKBB(2,3)= VALMAT(14) 694 : KKBB(3,1)= VALMAT(15) 695 : KKBB(3,2)= VALMAT(16) 696 : KKBB(3,3)= VALMAT(17) 697 : ENDIF 698 : c 698 : CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS) 699 : c 700 : c calcul des eps 2 701 : c 702 : IF (IREPS2.EQ.1) 702 : & CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,IGAU, 703 : & XDPGE,YDPGE,UDPGE,NHRM) 704 : c 705 : c contribution de epsi a msr0 706 : c 707 : TRACEP=XSTRS(1)+XSTRS(2)+XSTRS(3) 708 : DO 4873 K=1,IDECAP 709 : XSTRS(NSTRS-IDECAP+K)=CPBB(K)*TRACEP 710 : 4873 CONTINUE 711 : DO 4773 KA=1,LHOOK 712 : XWRK(KA)=XSTRS(KA) 713 : 4773 CONTINUE 714 : DO 7873 KA=1,LHOOK 715 : r_z = 0.D0 716 : DO KB=1,LHOOK 717 : r_z = r_z + DDHOOK(KA,KB)*XWRK(KB) 718 : ENDDO 719 : XSTRS(KA) = r_z 720 : 7873 CONTINUE 721 : c 722 : c calcul de l'effet de la pression 723 : c 723 : CALL SIGPOR(COBMA,UNSURM,XGENE,NSTN,XDDL,IFOUR,NSTRS, 724 : . XSTRS,LRB,LRN,LPP,MELE,I19,COBB,KKBB,IDECAP) 725 : c 726 : c remplissage du segment contenant les contraintes 727 : c 728 : MPTVAL=IVASTR 729 : DO 7073 ICOMP=1,NSTRS 730 : MELVAL=IVAL(ICOMP) 731 : IBMN=MIN(IB ,VELCHE(/2)) 732 : VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 733 : 7073 CONTINUE 734 : c 735 : 5173 CONTINUE 736 : c 737 : IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN 738 : INTERR(1)=IB 738 : CALL ERREUR(195) 739 : GOTO 9973 740 : ENDIF 741 : c 742 : 3173 CONTINUE 743 : c 743 : IF (I19.NE.0) CALL ERREUR(19) 744 : IF (IRTD.EQ.0) THEN 745 : MOTERR(1:8)=CMATE 746 : MOTERR(9:12)=NOMFR(MFR/2+1) 747 : INTERR(1)=IFOUR 747 : CALL ERREUR(81) 748 : ENDIF 749 : 9973 CONTINUE 750 : SEGSUP WRK2,WRK5 751 : IF (IREPS2.EQ.1) SEGSUP MTRACE 752 : GOTO 510 753 : 754 : c____________________________________________________________________ 755 : c 756 : c joints poreux 757 : c____________________________________________________________________ 758 : c 759 : 80 CONTINUE 760 : c Ces cas ne sont pas prevus actuellement ! 761 : IF ( IMAT.EQ.2 .OR. (IMAT.EQ.1.AND.MATE.NE.1) ) GOTO 99 762 : c 763 : c pour ces elements nbbb = nombre de noeuds 764 : c nbno = nombre de fonctions de forme 765 : c 766 : NBNO=IPORE 767 : NSTN=1 768 : NSTMU=2 769 : IF (IFOUR.EQ.2) NSTMU=3 770 : LRN=(NBNO-NBBB)*3/2 771 : LRB=LRE-LRN 772 : NFAC=(3*NBBB-NBNO)/2 773 : c 774 : SEGINI WRK2,WRK3,WRK5 775 : c 776 : I19 =0 777 : c 778 : DO 3080 IB=1,NBELEM 779 : c 780 : c on cherche d'abord les deplacements 781 : c 782 : MPTVAL=IVADEP 783 : IE=1 784 : DO 4180 IGAU=1,NFAC 785 : DO 4280 ICOMP=1,NDEP-1 786 : MELVAL=IVAL(ICOMP) 787 : IGMN=MIN(IGAU,VELCHE(/1)) 788 : IBMN=MIN(IB ,VELCHE(/2)) 789 : XDDL(IE)=VELCHE(IGMN,IBMN) 790 : IE=IE+1 791 : 4280 CONTINUE 792 : 4180 CONTINUE 793 : c 794 : c puis les pressions 795 : c 796 : MELVAL=IVAL(NDEP) 797 : DO 4080 IGAU=1,NBNN 798 : DO 4190 INSOM=1,NBSOM(IELE) 799 : IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GOTO 4191 800 : 4190 CONTINUE 801 : IF (IGAU.GT.NFAC) GOTO 4191 802 : GOTO 4080 803 : 4191 CONTINUE 804 : IBMN=MIN(IB ,VELCHE(/2)) 805 : IGMN=MIN(IGAU,VELCHE(/1)) 806 : XDDL(IE)=VELCHE(IGMN,IBMN) 807 : IE=IE+1 808 : 4080 CONTINUE 809 : c 810 : c on cherche les coordonnees des noeuds de l element ib 811 : c 811 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 812 : c 813 : c calcul des exes locaux et des coordonnees locales 814 : c 814 : CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS) 815 : c 816 : c boucle sur les points de gauss 817 : c 818 : ISDJC=0 818 : CALL ZERO(COBMA,LHOOK,1) 819 : c 820 : DO 5080 IGAU=1,NBPTEL 821 : c 821 : CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL, 822 : . SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,1) 823 : c 824 : IF (DJAC.EQ.0.D0) THEN 825 : INTERR(1)=IB 825 : CALL ERREUR(259) 826 : GOTO 9980 827 : ELSE IF (DJAC.LT.0.D0) THEN 828 : ISDJC=ISDJC+1 829 : ENDIF 830 : c 831 : MPTVAL=IVAMAT 832 : C*D IF (IMAT.EQ.2) THEN 833 : C*D GO TO 99 834 : C*D ELSE IF (IMAT.EQ.1) THEN 835 : DO 9080 IM=1,NMATT 836 : IF (IVAL(IM).NE.0) THEN 837 : MELVAL=IVAL(IM) 838 : IBMN=MIN(IB ,VELCHE(/2)) 839 : IGMN=MIN(IGAU,VELCHE(/1)) 840 : VALMAT(IM)=VELCHE(IGMN,IBMN) 841 : ELSE 842 : VALMAT(IM)=0.D0 843 : ENDIF 844 : 9080 CONTINUE 845 : C*D IF (MATE.EQ.1) THEN 846 : IF (IGAU.LE.NBGMAT.AND.(IB.LE.NELMAT.OR.NBGMAT.GT.1)) 846 : & CALL DOUO88(VALMAT,CMATE,IFOUR,LHOOK,DDHOOK,IRTD) 847 : C*D ELSE 848 : C*D GO TO 99 849 : C*D ENDIF 850 : C*D ENDIF 851 : c 851 : CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS) 852 : c 853 : c contribution de epsi a msr0 854 : c 855 : IF (IMAT.EQ.1) THEN 856 : COBMA(NSTMU)=VALMAT(3) 857 : XMOB=VALMAT(4) 858 : ENDIF 859 : XSTRS(NSTRS)=COBMA(NSTMU)*XSTRS(NSTMU) 860 : 861 : DO 4887 KA=1,LHOOK 862 : XWRK(KA)=XSTRS(KA) 863 : 4887 CONTINUE 864 : 865 : DO 4886 KA=1,LHOOK 866 : r_z = 0.D0 867 : DO KB=1,LHOOK 868 : r_z = r_z + DDHOOK(KA,KB)*XWRK(KB) 869 : ENDDO 870 : XSTRS(KA)= r_z 871 : 4886 CONTINUE 872 : c 873 : c calcul de l'effet de la pression 874 : c 875 : IF (XMOB.EQ.0.D0) THEN 876 : UNSURM=0.D0 877 : ELSE 878 : UNSURM=1.D0 / XMOB 879 : ENDIF 880 : c 880 : CALL SIGPOR(COBMA,UNSURM,XGENE,NSTN,XDDL,IFOUR,NSTRS, 881 : . XSTRS,LRB,LRN,LPP,MELE,I19,COBB,KKBB,IDECAP) 882 : c 883 : c remplissage du segment contenant les contraintes 884 : c 885 : MPTVAL=IVASTR 886 : DO 7080 ICOMP=1,NSTRS 887 : MELVAL=IVAL(ICOMP) 888 : IBMN=MIN(IB ,VELCHE(/2)) 889 : VELCHE(IGAU,IBMN)=XSTRS(ICOMP) 890 : 7080 CONTINUE 891 : c 892 : 5080 CONTINUE 893 : c 894 : IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN 895 : INTERR(1)=IB 895 : CALL ERREUR(195) 896 : GOTO 9980 897 : ENDIF 898 : c 899 : 3080 CONTINUE 900 : c 901 : IF (IRTD.EQ.0) THEN 902 : MOTERR(1:8)=CMATE 903 : MOTERR(9:12)=NOMFR(MFR/2+1) 904 : INTERR(1)=IFOUR 904 : CALL ERREUR(81) 905 : ENDIF 905 : IF (I19.NE.0) CALL ERREUR(19) 906 : 9980 CONTINUE 907 : SEGSUP,WRK2,WRK3,WRK5 908 : GOTO 510 909 : 910 : c____________________________________________________________________ 911 : 99 CONTINUE 912 : MOTERR(1:4)=NOMTP(MELE) 913 : MOTERR(9:12)='SIGM' 913 : CALL ERREUR(86) 914 : 915 : C- Fin du sous-programme 916 : 510 CONTINUE 917 : SEGSUP,MVELCH,WRK1 918 : 919 : RETURN 920 : END 921 : 922 : 923 : 924 :
© Cast3M 2003 - All rights reserved.
Disclaimer