1 : C FRIG2C SOURCE FANDEUR 11/10/07 21:15:21 7159 2 : 3 : SUBROUTINE FRIG2C (maifro,IPRIGI,IPCHJE,IPRIG2) 4 : 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8(A-H,O-Z) 7 : 8 : * Ce sous-programme calcule la raideur de frottement en 2D. 9 : * il a besoin pour cela du maillage de frottement et de la raideur 10 : * de contact (ou la raideur totale si c'est plus simple) 11 : 12 : -INC CCOPTIO 13 : -INC CCREEL 14 : 15 : -INC SMCHPOI 16 : -INC SMCOORD 17 : -INC SMELEME 18 : -INC SMRIGID 19 : 20 : * icpr lx du contact 21 : segment icpr(xcoor(/1)/idimp1) 22 : segment mwrk 23 : real*8 xdir(nbno),ydir(nbno) 24 : real*8 alpha(nbno),beta(nbno) 25 : integer in1(nbno),in2(nbno),in3(nbno),in4(nbno) 26 : endsegment 27 : 28 : character*4 modepl(4),moforc(4) 29 : data modepl /'UX ','UY ','UR ','UZ '/ 30 : data moforc /'FX ','FY ','FR ','FZ '/ 31 : * 32 : * test sur les options de calcul (faits dans frig) 33 : C* IF (IDIM .NE. 2) THEN 33 : C* CALL ERREUR(5) 34 : C* RETURN 35 : C* ENDIF 36 : C* IF (IFOMOD .NE. -1 .AND. IFOMOD .NE. 0) THEN 36 : C* CALL ERREUR(710) 37 : C* RETURN 38 : C* ENDIF 39 : * 40 : * quelques initialisations 41 : * 42 : IPRIG2 = 0 43 : idimp1 = IDIM + 1 44 : * indice servant a affecter le nom des composantes 45 : idepl = (ifomod + 1)*2 46 : * 47 : meleme = maifro 48 : SEGACT,meleme 49 : if (lisous(/1).ne.0) then 49 : call erreur(16) 50 : ELSE IF (itypel.NE.22) then 50 : CALL ERREUR(16) 51 : ELSE IF (num(/1).NE.5 .AND. num(/1).NE.6) THEN 51 : CALL ERREUR(16) 52 : ENDIF 53 : IF (IERR.NE.0) GOTO 9900 54 : 55 : * Lecture de la raideur de contact et extraction de tous les LX 56 : mwrk = 0 57 : segini icpr 58 : 59 : nbno=0 60 : nbel=0 61 : nbel3=0 62 : nbel4=0 63 : nbel5=0 64 : xcr=xpetit*1.E+10 65 : 66 : mrigid = IPRIGI 67 : SEGACT,mrigid 68 : 69 : DO 10 irig = 1, irigel(/2) 70 : IF (irigel(6,irig).EQ.0) GOTO 10 71 : ipt2 = irigel(1,irig) 72 : SEGACT,ipt2 73 : IF (ipt2.itypel.NE.22) GOTO 10 74 : ityp2 = ipt2.num(/1) 75 : nbel22 = ipt2.num(/2) 76 : nbel = nbel + nbel22 77 : c* Formulation forte/standard 78 : IF (ityp2.EQ.3) THEN 79 : nbel3 = nbel3 + nbel22 80 : ELSEIF (ityp2.EQ.4) THEN 81 : nbel4 = nbel4 + nbel22 82 : c* Formulation faible 83 : ELSEIF (ityp2.EQ.5) THEN 84 : nbel5 = nbel5 + nbel22 85 : ELSE 85 : CALL ERREUR(16) 86 : goto 10 87 : ENDIF 88 : i1 = (ipt2.num(2,1)-1)*idimp1 89 : x1 = xcoor(i1+1) 90 : y1 = xcoor(i1+2) 91 : do 101 iel = 1, nbel22 92 : nbno=nbno+1 93 : icpr(ipt2.num(1,iel))=nbno 94 : i2 = (ipt2.num(2,iel)-1)*idimp1 95 : xcr = MAX(xcr,ABS(x1-xcoor(i2+1)), 96 : & ABS(y1-xcoor(i2+2))) 97 : 101 continue 98 : 10 CONTINUE 99 : 100 : xcr=xcr*1.E-12 101 : * write (ioimp,*) ' xcr ',xcr 102 : 103 : segini mwrk 104 : 105 : DO 11 irig = 1, irigel(/2) 106 : if (irigel(6,irig).eq.0) goto 11 107 : ipt2 = irigel(1,irig) 108 : IF (ipt2.itypel.ne.22) goto 110 109 : ityp2 = ipt2.num(/1) 110 : nbel22 = ipt2.num(/2) 111 : xmatri=irigel(4,irig) 112 : segact xmatri 113 : c* Relation type 1 : noeud-noeud (forte/standard) 114 : if (ityp2.eq.3) then 115 : do 113 iel = 1, nbel22 116 : ip = icpr(ipt2.num(1,iel)) 117 : c* if (ip.eq.0) goto 113 118 : in1(ip)=ipt2.num(2,iel) 119 : in2(ip)=ipt2.num(3,iel) 120 : xdir(ip)=re(1,4,iel) 121 : ydir(ip)=re(1,5,iel) 122 : c* alpha(ip)=+1. 123 : c* write(ioimp,*) ' ip ',ip 124 : 113 continue 125 : c* Relation type 2 : noeud-segment (forte/standard) 126 : elseif (ityp2.eq.4) then 127 : do 114 iel = 1, nbel22 128 : ip = icpr(ipt2.num(1,iel)) 129 : c* if (ip.eq.0) goto 114 130 : in1(ip)=ipt2.num(2,iel) 131 : in2(ip)=ipt2.num(3,iel) 132 : in3(ip)=ipt2.num(4,iel) 133 : xdir(ip)=re(1,6,iel) 134 : ydir(ip)=re(1,7,iel) 135 : if (ABS(xdir(ip)).gt.0.1) then 136 : alpha(ip)=-re(1,2,iel)/xdir(ip) 137 : else 138 : alpha(ip)=-re(1,3,iel)/ydir(ip) 139 : endif 140 : beta(ip) = 1. - alpha(ip) 141 : c* write(ioimp,*) ' ip alpha ',ip,alpha(ip) 142 : 114 continue 143 : c* Relation type 0 : segment-segment (faible) 144 : elseif (ityp2.eq.5) then 145 : do 115 iel = 1, nbel22 146 : ip = icpr(ipt2.num(1,iel)) 147 : c* if (ip.eq.0) goto 115 148 : in1(ip)=ipt2.num(2,iel) 149 : in2(ip)=ipt2.num(3,iel) 150 : in3(ip)=ipt2.num(4,iel) 151 : in4(ip)=ipt2.num(5,iel) 152 : xdir(ip)=re(1,2,iel)+re(1,4,iel) 153 : ydir(ip)=re(1,3,iel)+re(1,5,iel) 154 : * Normale de norme = 1 155 : * xl=sqrt(xdir(ip)**2+ydir(ip)**2) 156 : * if (ABS(xdir(ip)/xl).gt.0.1) then 157 : if (ABS(xdir(ip)).gt.0.1) then 158 : alpha(ip)= re(1,2,iel)/xdir(ip) 159 : beta(ip) =-re(1,6,iel)/xdir(ip) 160 : else 161 : alpha(ip)= re(1,3,iel)/ydir(ip) 162 : beta(ip) =-re(1,7,iel)/ydir(ip) 163 : endif 164 : c* write(ioimp,*) ' ip alpha beta ',ip,alpha(ip),beta(ip) 165 : 115 continue 166 : endif 167 : segdes xmatri 168 : 110 continue 169 : segdes ipt2 170 : 11 CONTINUE 171 : segdes mrigid 172 : IF (ierr.NE.0) GOTO 9910 173 : 174 : * Traitement du champ de jeux (demi-frottement si jeu non nul) 175 : mchpoi = IPCHJE 176 : segact mchpoi 177 : do 15 isoupo = 1, ipchp(/1) 178 : msoupo = ipchp(isoupo) 179 : segact msoupo 180 : mpoval=ipoval 181 : segact mpoval 181 : if (vpocha(/2).ne.1) call erreur(16) 182 : ipt8=igeoc 183 : segact ipt8 183 : if (ipt8.num(/1).ne.vpocha(/2)) call erreur(16) 184 : do 16 i=1,vpocha(/1) 185 : if ((vpocha(i,1)).gt.xcr) then 186 : icpr(ipt8.num(1,i))=-ABS(icpr(ipt8.num(1,i))) 187 : * write(ioimp,*) ' frig2c frottement diminué ',ipt8.num(1,i), 188 : * & vpocha(i,1) 189 : else 190 : * write(ioimp,*) ' frig2c frottement gardé ',ipt8.num(1,i), 191 : * & vpocha(i,1) 192 : endif 193 : 16 continue 194 : segdes ipt8,mpoval,msoupo 195 : 15 continue 196 : segdes mchpoi 197 : IF (ierr.NE.0) GOTO 9910 198 : 199 : * creation des raideurs de frottement la ou il existe une raideur de 200 : * contact 201 : nrigel=0 202 : if (nbel3.ne.0) nrigel=nrigel+1 203 : if (nbel4.ne.0) nrigel=nrigel+1 204 : if (nbel5.ne.0) nrigel=nrigel+1 205 : * write(ioimp,*) ' nrigel ',nrigel,nbel3,nbel4,nbel5 206 : * 207 : segini mrigid 208 : * 209 : * Tout ce qu'on cree est du frottement et est symetrique 210 : DO irig = 1,nrigel 211 : irigel(6,irig) = 2 212 : irigel(7,irig) = 0 213 : coerig(irig)=1.D0 214 : enddo 215 : ichole=0 216 : imgeo1=0 217 : imgeo2=0 218 : isupeq=0 219 : iforig=ifomod 220 : mtymat='RIGIDITE' 221 : * 222 : irig=0 223 : * 224 : * Relation type 0 : segment-segment (faible) 225 : *------------------ 226 : if (nbel5.ne.0) then 227 : irig=irig+1 228 : nelri5=nbel5 229 : * 230 : * creation du meleme associe a la relation 231 : * 232 : nbelem=nelri5 233 : nbnn =5 234 : nbsous=0 235 : nbref =0 236 : segini ipt2 237 : ipt2.itypel=22 238 : irigel(1,irig)=ipt2 239 : * 240 : * creation du descriptif commun a toutes les raideurs 241 : * 242 : nligrp=9 243 : nligrd=nligrp 244 : segini descr 245 : lisinc(1)='LX ' 246 : lisdua(1)='FLX ' 247 : noelep(1)=1 248 : noeled(1)=1 249 : do 501 i=2,9,2 250 : lisinc(i )=modepl(1+idepl) 251 : lisinc(i+1)=modepl(2+idepl) 252 : lisdua(i )=moforc(1+idepl) 253 : lisdua(i+1)=moforc(2+idepl) 254 : noelep(i )=(i+2)/2 255 : noelep(i+1)=noelep(i) 256 : noeled(i )=noelep(i) 257 : noeled(i+1)=noelep(i) 258 : 501 continue 259 : segdes descr 260 : irigel(3,irig)=descr 261 : * 262 : * creation du segment xmatri 263 : * 264 : nelrig=nelri5 265 : segini xmatri 266 : irigel(4,irig)=xmatri 267 : * 268 : * on boucle sur le maillage de frottement en creant la raideur uniquement 269 : * la ou il y a une raideur de contact 270 : * 271 : nelri5=0 272 : * 273 : do 500 iel=1,num(/2) 274 : ip=icpr(num(1,iel)) 275 : if (ip.eq.0) goto 500 276 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 277 : if (ip.lt.0) then 278 : xcof=0.5d0 279 : else 280 : xcof=1.d0 281 : endif 282 : ip=ABS(ip) 283 : * 284 : * on range le maillage sous-jacent dans le meleme 285 : * 286 : ip4=in4(ip) 287 : if (ip4.eq.0) goto 500 288 : ip1 =in1(ip) 289 : ip2 =in2(ip) 290 : ip3 =in3(ip) 291 : nelri5=nelri5+1 292 : ipt2.num(1,nelri5)=num(6,iel) 293 : ipt2.num(2,nelri5)=ip1 294 : ipt2.num(3,nelri5)=ip2 295 : ipt2.num(4,nelri5)=ip3 296 : ipt2.num(5,nelri5)=ip4 297 : * 298 : xr = xcof*ydir(ip) 299 : yr = -xcof*xdir(ip) 300 : * 301 : * on remplit le xmatri 302 : * 303 : * lambda 304 : re(1,1,nelri5) = 0. 305 : * ip1 306 : re(2,1,nelri5) = alpha(ip)*xr 307 : re(3,1,nelri5) = alpha(ip)*yr 308 : * ip2 309 : re(4,1,nelri5) = (1.-alpha(ip))*xr 310 : re(5,1,nelri5) = (1.-alpha(ip))*yr 311 : * ip3 312 : re(6,1,nelri5) = -beta(ip)*xr 313 : re(7,1,nelri5) = -beta(ip)*yr 314 : * ip 315 : re(8,1,nelri5) = (beta(ip)-1.)*xr 316 : re(9,1,nelri5) = (beta(ip)-1.)*yr 317 : * on transpose 318 : do 510 ic=2,nligrd 319 : re(1,ic,nelri5)=re(ic,1,nelri5) 320 : 510 continue 321 : * 322 : 500 continue 323 : 324 : if (nelrig.ne.nelri5) then 325 : nelrig=nelri5 326 : segadj,xmatri 327 : nbelem=nelri5 328 : nbnn =5 329 : nbsous=0 330 : nbref =0 331 : segadj,ipt2 332 : endif 333 : segdes,xmatri,ipt2 334 : endif 335 : * 336 : * Relation type 2 : noeud-segment (forte/standard) 337 : * ----------------- 338 : if (nbel4.ne.0) then 339 : * 340 : irig=irig+1 341 : nelri4=nbel4 342 : * 343 : * creation du meleme associe a la relation 344 : * 345 : nbelem=nelri4 346 : nbnn =4 347 : nbsous=0 348 : nbref =0 349 : segini ipt2 350 : ipt2.itypel=22 351 : irigel(1,irig)=ipt2 352 : * 353 : * creation du descriptif commun a toutes les raideurs 354 : * 355 : nligrp=7 356 : nligrd=nligrp 357 : segini,descr 358 : lisinc(1)='LX ' 359 : lisdua(1)='FLX ' 360 : noelep(1)=1 361 : noeled(1)=1 362 : do 401 i=2,nligrp,2 363 : lisinc(i) =modepl(1+idepl) 364 : lisinc(i+1)=modepl(2+idepl) 365 : lisdua(i) =moforc(1+idepl) 366 : lisdua(i+1)=moforc(2+idepl) 367 : noelep(i) =(i+2)/2 368 : noelep(i+1)=noelep(i) 369 : noeled(i) =noelep(i) 370 : noeled(i+1)=noelep(i) 371 : 401 continue 372 : segdes,descr 373 : irigel(3,irig)=descr 374 : * 375 : * creation du segment xmatri 376 : * 377 : nelrig=nbel4 378 : segini xmatri 379 : irigel(4,irig)=xmatri 380 : * 381 : * on boucle sur le maillage de frottement en creant la raideur uniquement 382 : * la ou il y a une raideur de contact 383 : * 384 : nelri4 = 0 385 : * 386 : do 400 iel=1,num(/2) 387 : * 388 : ip=icpr(num(1,iel)) 389 : if (ip.eq.0) goto 400 390 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 391 : if (ip.lt.0) then 392 : xcof=0.5d0 393 : else 394 : xcof=1.d0 395 : endif 396 : ip=ABS(ip) 397 : * 398 : * on range le maillage sous-jacent dans le meleme 399 : * 400 : c* Tests a voir 401 : if (in4(ip).ne.0) goto 400 402 : ip3 =in3(ip) 403 : if (ip3.eq.0) goto 400 404 : ip1 =in1(ip) 405 : ip2 =in2(ip) 406 : nelri4=nelri4+1 407 : ipt2.num(1,nelri4)=num(5,iel) 408 : ipt2.num(2,nelri4)=ip1 409 : ipt2.num(3,nelri4)=ip2 410 : ipt2.num(4,nelri4)=ip3 411 : * 412 : xr = xcof*ydir(ip) 413 : yr = -xcof*xdir(ip) 414 : c* OF A voir !!!!!!!!!!!!!!!!!!! 415 : c* xr = -xcof*ydir(ip) 416 : c* yr = xcof*xdir(ip) 417 : * 418 : * on remplit le xmatri 419 : * 420 : * lambda 421 : re(1,1,nelri4) = 0. 422 : * ip1 423 : re(2,1,nelri4) = alpha(ip)*xr 424 : re(3,1,nelri4) = alpha(ip)*yr 425 : * ip2 426 : re(4,1,nelri4) = beta(ip)*xr 427 : re(5,1,nelri4) = beta(ip)*yr 428 : * ip 429 : re(6,1,nelri4) = -xr 430 : re(7,1,nelri4) = -yr 431 : * on transpose 432 : do 410 ic=2,nligrp 433 : re(1,ic,nelri4)=re(ic,1,nelri4) 434 : 410 continue 435 : * 436 : 400 continue 437 : 438 : if (nelrig.ne.nelri4) then 439 : nelrig=nelri4 440 : segadj,xmatri 441 : nbelem=nelri4 442 : nbsous=0 443 : nbref=0 444 : nbnn=4 445 : segadj,ipt2 446 : endif 447 : segdes xmatri,ipt2 448 : endif 449 : * 450 : * Relation type 1 : noeud-noeud (forte/standard) 451 : * 452 : if (nbel3.ne.0) then 453 : * 454 : irig=irig+1 455 : nelri3=nbel3 456 : * 457 : * creation du meleme associe a la relation 458 : * 459 : nbsous=0 460 : nbref=0 461 : nbnn=3 462 : nbelem=nelri3 463 : segini ipt2 464 : ipt2.itypel=22 465 : irigel(1,irig)=ipt2 466 : * 467 : * creation du descriptif commun a toutes les raideurs 468 : * 469 : nligrp=5 470 : nligrd=nligrp 471 : segini descr 472 : lisinc(1)='LX ' 473 : lisdua(1)='FLX ' 474 : noelep(1)=1 475 : noeled(1)=1 476 : do 301 i=2,nligrp,2 477 : lisinc(i) =modepl(1+idepl) 478 : lisinc(i+1)=modepl(2+idepl) 479 : lisdua(i) =moforc(1+idepl) 480 : lisdua(i+1)=moforc(2+idepl) 481 : noelep(i) =(i+2)/2 482 : noelep(i+1)=noelep(i) 483 : noeled(i) =noelep(i) 484 : noeled(i+1)=noelep(i) 485 : 301 continue 486 : segdes,descr 487 : irigel(3,irig)=descr 488 : * 489 : * creation du segment xmatri 490 : * 491 : nelrig=nelri3 492 : segini xmatri 493 : irigel(4,irig)=xmatri 494 : * 495 : * on boucle sur le maillage de frottement en creant la raideur uniquement 496 : * la ou il y a une raideur de contact 497 : * 498 : nelri3 = 0 499 : do 300 iel=1,num(/2) 500 : ip=icpr(num(1,iel)) 501 : if (ip.eq.0) goto 300 502 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 503 : if (ip.lt.0) then 504 : xcof=0.5d0 505 : else 506 : xcof=1.d0 507 : endif 508 : ip=ABS(ip) 509 : * 510 : * on range le maillage sous-jacent dans le meleme 511 : * 512 : if (in3(ip).ne.0) goto 300 513 : ip1 =in1(ip) 514 : ip2 =in2(ip) 515 : nelri3=nelri3+1 516 : ipt2.num(1,nelri3)=num(5,iel) 517 : ipt2.num(2,nelri3)=ip1 518 : ipt2.num(3,nelri3)=ip2 519 : * 520 : xr = xcof*ydir(ip) 521 : yr = -xcof*xdir(ip) 522 : * 523 : * on remplit le xmatri 524 : * 525 : * lambda 526 : re(1,1,nelri3)= 0. 527 : * ip1 528 : re(2,1,nelri3)= xr 529 : re(3,1,nelri3)= yr 530 : * ip 531 : re(4,1,nelri3)= -xr 532 : re(5,1,nelri3)= -yr 533 : * on transpose 534 : do 310 ic=2,nligrd 535 : re(1,ic,nelri3)=re(ic,1,nelri3) 536 : 310 continue 537 : * 538 : 300 continue 539 : * 540 : if (nelrig.ne.nelri3) then 541 : nelrig=nelri3 542 : segadj,xmatri 543 : nbelem=nelri3 544 : nbnn =3 545 : nbsous=0 546 : nbref =0 547 : segadj,ipt2 548 : endif 549 : segdes,xmatri,ipt2 550 : endif 551 : 552 : * segdes,mrigid 553 : IPRIG2 = mrigid 554 : 555 : 9910 CONTINUE 556 : SEGSUP,icpr,mwrk 557 : 9900 CONTINUE 558 : SEGDES,meleme 559 : 560 : RETURN 561 : END 562 : 563 :
© Cast3M 2003 - All rights reserved.
Disclaimer