1 : C FRIG3C SOURCE FANDEUR 11/10/07 21:15:23 7159 2 : 3 : SUBROUTINE FRIG3C (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 3D 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 SMMODEL 19 : -INC SMRIGID 20 : 21 : * icpr lx du contact 22 : segment icpr(xcoor(/1)/idimp1) 23 : segment mwrk 24 : real*8 xdir(nbno),ydir(nbno),zdir(nbno) 25 : real*8 alpha(nbno),beta(nbno),gamma(nbno) 26 : integer in1(nbno),in2(nbno),in3(nbno),in4(nbno) 27 : endsegment 28 : segment mwrkf 29 : real*8 delta(nbno) 30 : integer in5(nbno),in6(nbno) 31 : endsegment 32 : 33 : character*4 modepl(3),moforc(3) 34 : data modepl /'UX ','UY ','UZ '/ 35 : data moforc /'FX ','FY ','FZ '/ 36 : * 37 : * test sur les options de calcul (faits dans frig) 38 : C* IF (IDIM .NE. 3) THEN 38 : C* CALL ERREUR(5) 39 : C* RETURN 40 : C* ENDIF 41 : * 42 : * quelques initialisations 43 : * 44 : IPRIGI2 = 0 45 : * 46 : idimp1 = IDIM + 1 47 : * 48 : meleme = maifro 49 : segact,meleme 50 : if (lisous(/1).ne.0) then 50 : call erreur(16) 51 : GOTO 9900 52 : endif 52 : if (itypel.ne.22) call erreur(16) 52 : if (num(/1).ne.7 .and. num(/1).ne.9) call erreur(16) 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 : mwrkf = 0 58 : segini icpr 59 : 60 : nbno=0 61 : nbel=0 62 : nbel3=0 63 : nbel4=0 64 : nbel5=0 65 : nbel7=0 66 : xcr=xpetit*1.e+10 67 : 68 : mrigid = IPRIGI 69 : SEGACT,mrigid 70 : 71 : DO 10 irig = 1, irigel(/2) 72 : if (irigel(6,irig).eq.0) goto 10 73 : ipt2=irigel(1,irig) 74 : segact ipt2 75 : IF (ipt2.itypel.ne.22) goto 10 76 : ityp2 = ipt2.num(/1) 77 : nbel22 = ipt2.num(/2) 78 : nbel = nbel + nbel22 79 : c* Formulation forte/standard 80 : if (ityp2.EQ.3) THEN 81 : nbel3 = nbel3 + nbel22 82 : ELSEIF (ityp2.EQ.4) THEN 83 : nbel4 = nbel4 + nbel22 84 : ELSEIF (ityp2.EQ.5) THEN 85 : nbel5 = nbel5 + nbel22 86 : c* Formulation faible 87 : ELSEIF (ityp2.EQ.7) THEN 88 : nbel7 = nbel7 + nbel22 89 : ELSE 89 : CALL ERREUR(16) 90 : goto 10 91 : ENDIF 92 : i1 = (ipt2.num(2,1)-1)*idimp1 93 : x1 = xcoor(i1+1) 94 : y1 = xcoor(i1+2) 95 : z1 = xcoor(i1+3) 96 : do 101 iel = 1, nbel22 97 : nbno=nbno+1 98 : icpr(ipt2.num(1,iel)) = nbno 99 : i2 = (ipt2.num(2,iel)-1)*idimp1 100 : xcr = MAX(xcr,ABS(x1-xcoor(i2+1)), 101 : & ABS(y1-xcoor(i2+2)), 102 : & ABS(z1-xcoor(i2+3))) 103 : 101 continue 104 : 10 CONTINUE 105 : 106 : xcr=xcr*1.E-12 107 : * write (ioimp,*) ' xcr ',xcr 108 : 109 : segini,mwrk 110 : if (nbel7.ne.0) segini,mwrkf 111 : 112 : DO 11 irig=1,irigel(/2) 113 : if (irigel(6,irig).eq.0) goto 11 114 : ipt2=irigel(1,irig) 115 : if (ipt2.itypel.ne.22) goto 110 116 : ityp2 = ipt2.num(/1) 117 : nbel22 = ipt2.num(/2) 118 : xmatri=irigel(4,irig) 119 : segact xmatri 120 : c* Relation type 1 : noeud-noeud (forte/standard) 121 : if (ityp2.eq.3) then 122 : do 113 iel = 1, nbel22 123 : ip = icpr(ipt2.num(1,iel)) 124 : c* if (ip.eq.0) goto 113 125 : in1(ip)=ipt2.num(2,iel) 126 : in2(ip)=ipt2.num(3,iel) 127 : xdir(ip)= re(1,2,iel) 128 : ydir(ip)= re(1,3,iel) 129 : zdir(ip)= re(1,4,iel) 130 : c* xdir(ip)=-re(1,5,iel) 131 : c* ydir(ip)=-re(1,6,iel) 132 : c* zdir(ip)=-re(1,7,iel) 133 : c* alpha(ip) = 1. 134 : * write(ioimp,*) ' ip 1 ',ip,xdir(ip),ydir(ip),zdir(ip) 135 : 113 continue 136 : c* Relation type 2 : noeud-segment (forte/standard) 137 : elseif (ityp2.eq.4) then 138 : do 114 iel = 1, nbel22 139 : ip = icpr(ipt2.num(1,iel)) 140 : c* if (ip.eq.0) goto 114 141 : in1(ip)=ipt2.num(2,iel) 142 : in2(ip)=ipt2.num(3,iel) 143 : in3(ip)=ipt2.num(4,iel) 144 : xdir(ip)=-re(1, 8,iel) 145 : ydir(ip)=-re(1, 9,iel) 146 : zdir(ip)=-re(1,10,iel) 147 : if (abs(xdir(ip)).gt.0.1) then 148 : alpha(ip)= re(1,2,iel)/xdir(ip) 149 : beta(ip) = re(1,5,iel)/xdir(ip) 150 : elseif (abs(ydir(ip)).gt.0.1) then 151 : alpha(ip)= re(1,3,iel)/ydir(ip) 152 : beta(ip) = re(1,6,iel)/ydir(ip) 153 : else 154 : alpha(ip)= re(1,4,iel)/zdir(ip) 155 : beta(ip) = re(1,7,iel)/zdir(ip) 156 : endif 157 : c* beta(ip) = 1. - alpha(ip) 158 : * write(ioimp,*) ' ip 2 ',ip,xdir(ip),ydir(ip),zdir(ip), 159 : * > alpha(ip),beta(ip) 160 : 114 continue 161 : c* Relation type 3 : noeud-triangle (forte/standard) 162 : elseif (ityp2.eq.5) then 163 : do 115 iel = 1, nbel22 164 : ip = icpr(ipt2.num(1,iel)) 165 : c* if (ip.eq.0) goto 115 166 : in1(ip)=ipt2.num(2,iel) 167 : in2(ip)=ipt2.num(3,iel) 168 : in3(ip)=ipt2.num(4,iel) 169 : in4(ip)=ipt2.num(5,iel) 170 : xdir(ip)=-re(1,11,iel) 171 : ydir(ip)=-re(1,12,iel) 172 : zdir(ip)=-re(1,13,iel) 173 : if (abs(xdir(ip)).gt.0.1) then 174 : alpha(ip)= re(1,2,iel)/xdir(ip) 175 : beta(ip) = re(1,5,iel)/xdir(ip) 176 : gamma(ip)= re(1,8,iel)/xdir(ip) 177 : elseif (abs(ydir(ip)).gt.0.1) then 178 : alpha(ip)= re(1,3,iel)/ydir(ip) 179 : beta(ip) = re(1,6,iel)/ydir(ip) 180 : gamma(ip)= re(1,9,iel)/ydir(ip) 181 : else 182 : alpha(ip)= re(1, 4,iel)/zdir(ip) 183 : beta(ip) = re(1, 7,iel)/zdir(ip) 184 : gamma(ip)= re(1,10,iel)/zdir(ip) 185 : endif 186 : c* gamma(ip)= 1. - beta(ip) - gamma(ip) 187 : * write (ioimp,*) ' ip 3 ',ip,xdir(ip),ydir(ip),zdir(ip), 188 : * > alpha(ip),beta(ip),gamma(ip) 189 : 115 continue 190 : c* Relation type 0 : triangle-triangle (faible) 191 : elseif (ityp2.eq.7) then 192 : do 117 iel = 1, nbel22 193 : ip = icpr(ipt2.num(1,iel)) 194 : c* if (ip.eq.0) goto 117 195 : in1(ip)=ipt2.num(2,iel) 196 : in2(ip)=ipt2.num(3,iel) 197 : in3(ip)=ipt2.num(4,iel) 198 : in4(ip)=ipt2.num(5,iel) 199 : in5(ip)=ipt2.num(6,iel) 200 : in6(ip)=ipt2.num(7,iel) 201 : xdir(ip)=re(1, 2,iel)+re(1, 5,iel)+re(1, 8,iel) 202 : ydir(ip)=re(1, 3,iel)+re(1, 6,iel)+re(1, 9,iel) 203 : zdir(ip)=re(1, 4,iel)+re(1, 7,iel)+re(1,10,iel) 204 : if (abs(xdir(ip)).gt.0.1) then 205 : alpha(ip)= re(1, 2,iel)/xdir(ip) 206 : beta(ip) = re(1, 5,iel)/xdir(ip) 207 : gamma(ip)=-re(1,11,iel)/xdir(ip) 208 : delta(ip)=-re(1,14,iel)/xdir(ip) 209 : elseif (abs(ydir(ip)).gt.0.1) then 210 : alpha(ip)= re(1, 3,iel)/ydir(ip) 211 : beta(ip) = re(1, 6,iel)/ydir(ip) 212 : gamma(ip)=-re(1,12,iel)/ydir(ip) 213 : delta(ip)=-re(1,15,iel)/ydir(ip) 214 : else 215 : alpha(ip)= re(1, 4,iel)/zdir(ip) 216 : beta(ip) = re(1, 7,iel)/zdir(ip) 217 : gamma(ip)=-re(1,13,iel)/zdir(ip) 218 : delta(ip)=-re(1,16,iel)/zdir(ip) 219 : endif 220 : * write (ioimp,*) ' ip 0 ',ip,xdir(ip),ydir(ip),zdir(ip), 221 : * > alpha(ip),beta(ip),gamma(ip),delta(ip) 222 : 117 continue 223 : endif 224 : segdes xmatri 225 : 110 continue 226 : segdes ipt2 227 : 11 CONTINUE 228 : segdes mrigid 229 : IF (ierr.NE.0) GOTO 9910 230 : 231 : * lecture du champ de jeux (pas de frottement si jeu non nul) 232 : mchpoi = IPCHJE 233 : segact mchpoi 234 : do 15 isoupo=1,ipchp(/1) 235 : msoupo=ipchp(isoupo) 236 : segact msoupo 237 : mpoval=ipoval 238 : segact mpoval 238 : if (vpocha(/2).ne.1) call erreur(16) 239 : ipt8=igeoc 240 : segact ipt8 240 : if (ipt8.num(/1).ne.vpocha(/2)) call erreur(16) 241 : do 16 i=1,vpocha(/1) 242 : if ((vpocha(i,1)).gt.xcr) then 243 : icpr(ipt8.num(1,i))=-abs(icpr(ipt8.num(1,i))) 244 : * write(ioimp,*) ' frig3c frottement diminué ',ipt8.num(1,i), 245 : * & vpocha(i,1) 246 : else 247 : * write(ioimp,*) ' frig3c frottement gardé ',ipt8.num(1,i), 248 : * & vpocha(i,1) 249 : endif 250 : 16 continue 251 : segdes ipt8,mpoval,msoupo 252 : 15 continue 253 : segdes mchpoi 254 : IF (ierr.NE.0) GOTO 9910 255 : * 256 : * creation des raideurs de frottement la ou il existe une raideur de 257 : * contact 258 : * 259 : nrigel=0 260 : if (nbel3.ne.0) nrigel=nrigel+1 261 : if (nbel4.ne.0) nrigel=nrigel+1 262 : if (nbel5.ne.0) nrigel=nrigel+1 263 : if (nbel7.ne.0) nrigel=nrigel+1 264 : * write(ioimp,*) ' nrigel ',nrigel,nbel3,nbel4,nbel5,nbel7 265 : * 266 : segini mrigid 267 : * 268 : ichole=0 269 : imgeo1=0 270 : imgeo2=0 271 : isupeq=0 272 : iforig=ifomod 273 : mtymat='RIGIDITE' 274 : * tout ce qu'on cree est du frottement et est symetrique 275 : do irig = 1, nrigel 276 : irigel(6,irig) = 2 277 : irigel(7,irig) = 0 278 : coerig(irig) = 1.D0 279 : enddo 280 : 281 : irig=0 282 : * 283 : * Relation type 3 : noeud-triangle (forte/standard) 284 : * ----------------- 285 : if (nbel5.ne.0) then 286 : * 287 : irig=irig+1 288 : nelri5=nbel5*2 289 : * 290 : * creation du meleme associe a la relation 291 : * 292 : nbelem=nelri5 293 : nbnn =5 294 : nbsous=0 295 : nbref =0 296 : segini ipt2 297 : ipt2.itypel=22 298 : irigel(1,irig)=ipt2 299 : * 300 : * creation du descriptif commun a toutes les raideurs 301 : * 302 : nligrp=13 303 : nligrd=nligrp 304 : segini descr 305 : lisinc(1)='LX ' 306 : lisdua(1)='FLX ' 307 : noelep(1)=1 308 : noeled(1)=1 309 : do 501 i=2,nligrp,3 310 : lisinc(i )=modepl(1) 311 : lisinc(i+1)=modepl(2) 312 : lisinc(i+2)=modepl(3) 313 : lisdua(i )=moforc(1) 314 : lisdua(i+1)=moforc(2) 315 : lisdua(i+2)=moforc(3) 316 : noelep(i )=(i+4)/3 317 : noelep(i+1)=noelep(i) 318 : noelep(i+2)=noelep(i) 319 : noeled(i )=noelep(i) 320 : noeled(i+1)=noelep(i) 321 : noeled(i+2)=noelep(i) 322 : 501 continue 323 : segdes,descr 324 : irigel(3,irig)=descr 325 : * 326 : * creation du segment xmatri 327 : * 328 : nelrig=nelri5 329 : segini xmatri 330 : irigel(4,irig)=xmatri 331 : * 332 : * on boucle sur le maillage de frottement en creant les raideurs uniquement 333 : * la ou il y a une raideur de contact 334 : * 335 : nelri5=0 336 : * 337 : do 500 iel = 1, num(/2) 338 : * 339 : ip=icpr(num(1,iel)) 340 : if (ip.eq.0) goto 500 341 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 342 : if (ip.lt.0) then 343 : xcof=0.5d0 344 : else 345 : xcof=1.d0 346 : endif 347 : ip=ABS(ip) 348 : * 349 : * on range le maillage sous-jacent dans le meleme 350 : * 351 : ip4 =in4(ip) 352 : c* test ci-dessous ? 353 : if (ip4.eq.0) goto 500 354 : ip3 =in3(ip) 355 : ip2 =in2(ip) 356 : ip1 =in1(ip) 357 : nelri5=nelri5+1 358 : ipt2.num(1,nelri5)=num(6,iel) 359 : ipt2.num(2,nelri5)=ip1 360 : ipt2.num(3,nelri5)=ip2 361 : ipt2.num(4,nelri5)=ip3 362 : ipt2.num(5,nelri5)=ip4 363 : melri5=nelri5 364 : nelri5=nelri5+1 365 : ipt2.num(1,nelri5)=num(7,iel) 366 : ipt2.num(2,nelri5)=ip1 367 : ipt2.num(3,nelri5)=ip2 368 : ipt2.num(4,nelri5)=ip3 369 : ipt2.num(5,nelri5)=ip4 370 : * 371 : * on remplit les xmatri 372 : * 373 : * directions orthogonales à xdir ydir zdir 374 : * 375 : xr=-2.*zdir(ip) 376 : yr=-3.*zdir(ip) 377 : zr= 2.*xdir(ip)+3.*ydir(ip) 378 : if (ABS(xr).le.0.1.AND.ABS(yr).le.0.1.AND.ABS(zr).LE.0.1) then 379 : xr=-2.*ydir(ip) 380 : yr= 2.*xdir(ip)+3.*zdir(ip) 381 : zr=-3.*ydir(ip) 382 : endif 383 : sr=SQRT(xr*xr+yr*yr+zr*zr) 384 : xr=xr/sr 385 : yr=yr/sr 386 : zr=zr/sr 387 : xs=ydir(ip)*zr-zdir(ip)*yr 388 : ys=zdir(ip)*xr-xdir(ip)*zr 389 : zs=xdir(ip)*yr-ydir(ip)*xr 390 : ss=SQRT(xs*xs+ys*ys+zs*zs) 391 : xs=xs/ss 392 : ys=ys/ss 393 : zs=zs/ss 394 : * 395 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 396 : * corrections sur les directions (xyz)r et (xyz)s 397 : alp_d = alpha(ip)*xcof 398 : bet_d = beta(ip)*xcof 399 : gam_d = gamma(ip)*xcof 400 : * 401 : * lambda 402 : re(1,1,melri5) = 0. 403 : re(1,1,nelri5) = 0. 404 : * ip1 405 : re(2,1,melri5) = alp_d*xr 406 : re(3,1,melri5) = alp_d*yr 407 : re(4,1,melri5) = alp_d*zr 408 : re(2,1,nelri5) = alp_d*xs 409 : re(3,1,nelri5) = alp_d*ys 410 : re(4,1,nelri5) = alp_d*zs 411 : * ip2 412 : re(5,1,melri5) = bet_d*xr 413 : re(6,1,melri5) = bet_d*yr 414 : re(7,1,melri5) = bet_d*zr 415 : re(5,1,nelri5) = bet_d*xs 416 : re(6,1,nelri5) = bet_d*ys 417 : re(7,1,nelri5) = bet_d*zs 418 : * ip3 419 : re( 8,1,melri5) = gam_d*xr 420 : re( 9,1,melri5) = gam_d*yr 421 : re(10,1,melri5) = gam_d*zr 422 : re( 8,1,nelri5) = gam_d*xs 423 : re( 9,1,nelri5) = gam_d*ys 424 : re(10,1,nelri5) = gam_d*zs 425 : * jp = ip4 426 : re(11,1,melri5) = -xr*xcof 427 : re(12,1,melri5) = -yr*xcof 428 : re(13,1,melri5) = -zr*xcof 429 : re(11,1,nelri5) = -xs*xcof 430 : re(12,1,nelri5) = -ys*xcof 431 : re(13,1,nelri5) = -zs*xcof 432 : * on transpose 433 : do 510 ic=2,nligrp 434 : re(1,ic,melri5)=re(ic,1,melri5) 435 : re(1,ic,nelri5)=re(ic,1,nelri5) 436 : 510 continue 437 : * 438 : 500 continue 439 : * 440 : if (nelrig.ne.nelri5) then 441 : nelrig=nelri5 442 : segadj,xmatri 443 : nbelem=nelri5 444 : nbnn =5 445 : nbsous=0 446 : nbref =0 447 : segadj,ipt2 448 : endif 449 : * 450 : segdes,xmatri,ipt2 451 : endif 452 : * 453 : * Relation type 2 : noeud-segment (forte/standard) 454 : * ----------------- 455 : if (nbel4.ne.0) then 456 : * 457 : irig=irig+1 458 : nelri4 = nbel4*2 459 : 460 : * creation du meleme associe a la relation 461 : * 462 : nbelem=nelri4 463 : nbnn =4 464 : nbsous=0 465 : nbref =0 466 : segini ipt2 467 : ipt2.itypel=22 468 : irigel(1,irig)=ipt2 469 : * 470 : * creation du descriptif commun a toutes les raideurs 471 : * 472 : nligrp=10 473 : nligrd=nligrp 474 : segini,descr 475 : lisinc(1)='LX ' 476 : lisdua(1)='FLX ' 477 : noelep(1)=1 478 : noeled(1)=1 479 : do 401 i=2,nligrp,3 480 : lisinc(i )=modepl(1) 481 : lisinc(i+1)=modepl(2) 482 : lisinc(i+2)=modepl(3) 483 : lisdua(i )=moforc(1) 484 : lisdua(i+1)=moforc(2) 485 : lisdua(i+2)=moforc(3) 486 : noelep(i )=(i+4)/3 487 : noelep(i+1)=noelep(i) 488 : noelep(i+2)=noelep(i) 489 : noeled(i )=noelep(i) 490 : noeled(i+1)=noelep(i) 491 : noeled(i+2)=noelep(i) 492 : 401 continue 493 : segdes,descr 494 : irigel(3,irig)=descr 495 : * 496 : * creation du segment xmatri 497 : * 498 : nelrig = nelri4 499 : segini xmatri 500 : irigel(4,irig)=xmatri 501 : * 502 : * on boucle sur le maillage de frottement en creant les raideurs uniquement 503 : * la ou il y a une raideur de contact 504 : * 505 : nelri4=0 506 : * 507 : do 400 iel=1,num(/2) 508 : ip=icpr(num(1,iel)) 509 : if (ip.eq.0) goto 400 510 : * coefficient multiplicateur selon signe de ip ("frottement diminue") 511 : if (ip.lt.0) then 512 : xcof=0.5d0 513 : else 514 : xcof=1.d0 515 : endif 516 : ip=ABS(ip) 517 : * 518 : * on range le maillage sous-jacent dans le meleme 519 : * 520 : c* Voir les tests ci-dessous ? 521 : if (in4(ip).ne.0) goto 400 522 : ip3 =in3(ip) 523 : if (ip3.eq.0) goto 400 524 : ip1 =in1(ip) 525 : ip2 =in2(ip) 526 : nelri4=nelri4+1 527 : ipt2.num(1,nelri4)=num(6,iel) 528 : ipt2.num(2,nelri4)=ip1 529 : ipt2.num(3,nelri4)=ip2 530 : ipt2.num(4,nelri4)=ip3 531 : melri4=nelri4 532 : nelri4=nelri4+1 533 : ipt2.num(1,nelri4)=num(7,iel) 534 : ipt2.num(2,nelri4)=ip1 535 : ipt2.num(3,nelri4)=ip2 536 : ipt2.num(4,nelri4)=ip3 537 : * 538 : * on remplit les xmatri 539 : * 540 : * directions orthogonales à xdir ydir zdir 541 : * 542 : xr=-2.*zdir(ip) 543 : yr=-3.*zdir(ip) 544 : zr=2.*xdir(ip)+3.*ydir(ip) 545 : if (abs(xr).le.0.1.and.abs(yr).le.0.1.and.abs(zr).le.0.1) then 546 : xr=-2.*ydir(ip) 547 : yr=2.*xdir(ip)+3.*zdir(ip) 548 : zr=-3.*ydir(ip) 549 : endif 550 : sr=sqrt(xr*xr+yr*yr+zr*zr) 551 : xr=xr/sr 552 : yr=yr/sr 553 : zr=zr/sr 554 : xs=ydir(ip)*zr-zdir(ip)*yr 555 : ys=zdir(ip)*xr-xdir(ip)*zr 556 : zs=xdir(ip)*yr-ydir(ip)*xr 557 : ss=sqrt(xs*xs+ys*ys+zs*zs) 558 : xs=xs/ss 559 : ys=ys/ss 560 : zs=zs/ss 561 : * 562 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 563 : * corrections sur les directions (xyz)r et (xyz)s 564 : alp_d = alpha(ip)*xcof 565 : bet_d = beta(ip)*xcof 566 : * lambda 567 : re(1,1,melri4) = 0. 568 : re(1,1,nelri4) = 0. 569 : * ip1 570 : re(2,1,melri4) = alp_d*xr 571 : re(3,1,melri4) = alp_d*yr 572 : re(4,1,melri4) = alp_d*zr 573 : re(2,1,nelri4) = alp_d*xs 574 : re(3,1,nelri4) = alp_d*ys 575 : re(4,1,nelri4) = alp_d*zs 576 : * ip2 577 : re(5,1,melri4) = bet_d*xr 578 : re(6,1,melri4) = bet_d*yr 579 : re(7,1,melri4) = bet_d*zr 580 : re(5,1,nelri4) = bet_d*xs 581 : re(6,1,nelri4) = bet_d*ys 582 : re(7,1,nelri4) = bet_d*zs 583 : * jp = ip3 584 : re( 8,1,melri4) = -xr*xcof 585 : re( 9,1,melri4) = -yr*xcof 586 : re(10,1,melri4) = -zr*xcof 587 : re( 8,1,nelri4) = -xs*xcof 588 : re( 9,1,nelri4) = -ys*xcof 589 : re(10,1,nelri4) = -zs*xcof 590 : * 591 : * on transpose 592 : do 410 ic=2,nligrp 593 : re(1,ic,melri4)=re(ic,1,melri4) 594 : re(1,ic,nelri4)=re(ic,1,nelri4) 595 : 410 continue 596 : 400 continue 597 : * 598 : if (nelrig.ne.nelri4) then 599 : nelrig=nelri4 600 : segadj,xmatri 601 : nbelem=nelri4 602 : nbnn =4 603 : nbsous=0 604 : nbref =0 605 : segadj,ipt2 606 : endif 607 : * 608 : segdes xmatri,ipt2,descr 609 : endif 610 : * 611 : * Relation type 1 : noeud-noeud (forte/standard) 612 : * ----------------- 613 : if (nbel3.ne.0) then 614 : * 615 : irig=irig+1 616 : nelri3 = nbel3*2 617 : * 618 : * creation du meleme associe a la relation 619 : * 620 : nbelem = nelri3 621 : nbnn = 3 622 : nbsous = 0 623 : nbref = 0 624 : segini ipt2 625 : ipt2.itypel=22 626 : irigel(1,irig)=ipt2 627 : * 628 : * creation du descriptif commun a toutes les raideurs 629 : * 630 : nligrp=7 631 : nligrd=nligrp 632 : segini,descr 633 : lisinc(1)='LX ' 634 : lisdua(1)='FLX ' 635 : noelep(1)=1 636 : noeled(1)=1 637 : do 301 i=2,nligrp,3 638 : lisinc(i )=modepl(1) 639 : lisinc(i+1)=modepl(2) 640 : lisinc(i+2)=modepl(3) 641 : lisdua(i )=moforc(1) 642 : lisdua(i+1)=moforc(2) 643 : lisdua(i+2)=moforc(3) 644 : noelep(i )=(i+4)/3 645 : noelep(i+1)=noelep(i) 646 : noelep(i+2)=noelep(i) 647 : noeled(i )=noelep(i) 648 : noeled(i+1)=noelep(i) 649 : noeled(i+2)=noelep(i) 650 : 301 continue 651 : segdes,descr 652 : irigel(3,irig)=descr 653 : * 654 : * creation du segment xmatri 655 : * 656 : nelrig = nelri3 657 : segini xmatri 658 : irigel(4,irig)=xmatri 659 : * 660 : * on boucle sur le maillage de frottement en creant les raideurs uniquement 661 : * la ou il y a une raideur de contact 662 : * 663 : nelri3=0 664 : * 665 : do 300 iel=1,num(/2) 666 : ip=icpr(num(1,iel)) 667 : if (ip.eq.0) goto 300 668 : * coefficient multiplicateur selon signe de ip ("frottement diminue") 669 : if (ip.lt.0) then 670 : xcof=0.5d0 671 : else 672 : xcof=1.d0 673 : endif 674 : ip=ABS(ip) 675 : * 676 : * on range le maillage sous-jacent dans le meleme 677 : * 678 : c* Voir les tests ci-dessous ? 679 : if (in3(ip).ne.0) goto 300 680 : if (in4(ip).ne.0) goto 300 681 : ip1 =in1(ip) 682 : ip2 =in2(ip) 683 : nelri3=nelri3+1 684 : ipt2.num(1,nelri3)=num(6,iel) 685 : ipt2.num(2,nelri3)=ip1 686 : ipt2.num(3,nelri3)=ip2 687 : melri3=nelri3 688 : nelri3=nelri3+1 689 : ipt2.num(1,nelri3)=num(7,iel) 690 : ipt2.num(2,nelri3)=ip1 691 : ipt2.num(3,nelri3)=ip2 692 : * 693 : * on remplit les xmatri 694 : * 695 : * directions orthogonales à xdir ydir zdir 696 : * 697 : xr=-2.*zdir(ip) 698 : yr=-3.*zdir(ip) 699 : zr=2.*xdir(ip)+3.*ydir(ip) 700 : if (abs(xr).le.0.1.and.abs(yr).le.0.1.and.abs(zr).le.0.1) then 701 : xr=-2.*ydir(ip) 702 : yr=2.*xdir(ip)+3.*zdir(ip) 703 : zr=-3.*ydir(ip) 704 : endif 705 : sr=sqrt(xr*xr+yr*yr+zr*zr) 706 : xr=xr/sr 707 : yr=yr/sr 708 : zr=zr/sr 709 : xs=ydir(ip)*zr-zdir(ip)*yr 710 : ys=zdir(ip)*xr-xdir(ip)*zr 711 : zs=xdir(ip)*yr-ydir(ip)*xr 712 : ss=sqrt(xs*xs+ys*ys+zs*zs) 713 : xs=xs/ss 714 : ys=ys/ss 715 : zs=zs/ss 716 : * 717 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 718 : * corrections sur les directions (xyz)r et (xyz)s 719 : xr = xr*xcof 720 : yr = yr*xcof 721 : zr = zr*xcof 722 : xs = xs*xcof 723 : ys = ys*xcof 724 : zs = zs*xcof 725 : * lambda 726 : re(1,1,melri3)= 0. 727 : re(1,1,nelri3)= 0. 728 : * ip1 729 : re(2,1,melri3)= xr 730 : re(3,1,melri3)= yr 731 : re(4,1,melri3)= zr 732 : re(2,1,nelri3)= xs 733 : re(3,1,nelri3)= ys 734 : re(4,1,nelri3)= zs 735 : * jp = ip2 736 : re(5,1,melri3)= -xr 737 : re(6,1,melri3)= -yr 738 : re(7,1,melri3)= -zr 739 : re(5,1,nelri3)= -xs 740 : re(6,1,nelri3)= -ys 741 : re(7,1,nelri3)= -zs 742 : * on transpose 743 : do 310 ic=2,nligrp 744 : re(1,ic,melri3)=re(ic,1,melri3) 745 : re(1,ic,nelri3)=re(ic,1,nelri3) 746 : 310 continue 747 : * 748 : 300 continue 749 : * 750 : if (nelrig.ne.nelri3) then 751 : nelrig=nelri3 752 : segadj,xmatri 753 : nbelem=nelri3 754 : nbnn =3 755 : nbsous=0 756 : nbref =0 757 : segadj,ipt2 758 : endif 759 : * 760 : segdes xmatri,ipt2 761 : endif 762 : * 763 : * Relation type 0 : triangle-triangle (faible) 764 : * ----------------- 765 : if (nbel7.ne.0) then 766 : * 767 : irig=irig+1 768 : nelri7=nbel7*2 769 : * 770 : * creation du meleme associe a la relation 771 : * 772 : nbelem=nelri7 773 : nbnn =7 774 : nbsous=0 775 : nbref =0 776 : segini ipt2 777 : ipt2.itypel=22 778 : irigel(1,irig)=ipt2 779 : * 780 : * creation du descriptif commun a toutes les raideurs 781 : * nligrp=idim*(nbnn-1)+1 782 : * 783 : nligrp=19 784 : nligrd=nligrp 785 : segini,descr 786 : lisinc(1)='LX ' 787 : lisdua(1)='FLX ' 788 : noelep(1)=1 789 : noeled(1)=1 790 : do 701 i=2,nligrp,3 791 : lisinc(i )=modepl(1) 792 : lisinc(i+1)=modepl(2) 793 : lisinc(i+2)=modepl(3) 794 : lisdua(i )=moforc(1) 795 : lisdua(i+1)=moforc(2) 796 : lisdua(i+2)=moforc(3) 797 : noelep(i )=(i+4)/3 798 : noelep(i+1)=noelep(i) 799 : noelep(i+2)=noelep(i) 800 : noeled(i )=noelep(i) 801 : noeled(i+1)=noelep(i) 802 : noeled(i+2)=noelep(i) 803 : 701 continue 804 : segdes,descr 805 : irigel(3,irig)=descr 806 : * 807 : * creation du segment xmatri 808 : * 809 : nelrig=nelri7 810 : segini xmatri 811 : irigel(4,irig)=xmatri 812 : * 813 : * on boucle sur le maillage de frottement en creant les raideurs uniquement 814 : * la ou il y a une raideur de contact 815 : * 816 : nelri7=0 817 : do 700 iel = 1, num(/2) 818 : ip=icpr(num(1,iel)) 819 : if (ip.eq.0) goto 700 820 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 821 : if (ip.lt.0) then 822 : xcof=0.5d0 823 : else 824 : xcof=1.d0 825 : endif 826 : ip=ABS(ip) 827 : * 828 : * on range le maillage sous-jacent dans le meleme 829 : * 830 : c* Voir les tests ci-dessous ? 831 : ip4 = in4(ip) 832 : if (ip4.eq.0) goto 700 833 : ip5 = in5(ip) 834 : if (ip5.eq.0) goto 700 835 : ip6 = in6(ip) 836 : if (ip6.eq.0) goto 700 837 : ip1 = in1(ip) 838 : ip2 = in2(ip) 839 : ip3 = in3(ip) 840 : nelri7=nelri7+1 841 : ipt2.num(1,nelri7)=num(8,iel) 842 : ipt2.num(2,nelri7)=ip1 843 : ipt2.num(3,nelri7)=ip2 844 : ipt2.num(4,nelri7)=ip3 845 : ipt2.num(5,nelri7)=ip4 846 : ipt2.num(6,nelri7)=ip5 847 : ipt2.num(7,nelri7)=ip6 848 : melri7=nelri7 849 : nelri7=nelri7+1 850 : ipt2.num(1,nelri7)=num(9,iel) 851 : ipt2.num(2,nelri7)=ip1 852 : ipt2.num(3,nelri7)=ip2 853 : ipt2.num(4,nelri7)=ip3 854 : ipt2.num(5,nelri7)=ip4 855 : ipt2.num(6,nelri7)=ip5 856 : ipt2.num(7,nelri7)=ip6 857 : * 858 : * on remplit les xmatri 859 : * 860 : * directions orthogonales à xdir ydir zdir 861 : * 862 : xr=-2.*zdir(ip) 863 : yr=-3.*zdir(ip) 864 : zr= 2.*xdir(ip)+3.*ydir(ip) 865 : if (ABS(xr).le.0.1.AND.ABS(yr).le.0.1.AND.ABS(zr).LE.0.1) then 866 : xr=-2.*ydir(ip) 867 : yr= 2.*xdir(ip)+3.*zdir(ip) 868 : zr=-3.*ydir(ip) 869 : endif 870 : sr=SQRT(xr*xr+yr*yr+zr*zr) 871 : xr=xr/sr 872 : yr=yr/sr 873 : zr=zr/sr 874 : xs=ydir(ip)*zr-zdir(ip)*yr 875 : ys=zdir(ip)*xr-xdir(ip)*zr 876 : zs=xdir(ip)*yr-ydir(ip)*xr 877 : ss=SQRT(xs*xs+ys*ys+zs*zs) 878 : xs=xs/ss 879 : ys=ys/ss 880 : zs=zs/ss 881 : * 882 : * coefficient multiplicateur selon signe de icpr ("frottement diminue") 883 : * corrections sur les directions (xyz)r et (xyz)s 884 : alp_d = alpha(ip)*xcof 885 : bet_d = beta(ip)*xcof 886 : uab_d = xcof - alp_d - bet_d 887 : gam_d = -gamma(ip)*xcof 888 : del_d = -delta(ip)*xcof 889 : ugd_d = -(xcof + gam_d + del_d) 890 : * 891 : * lambda 892 : re( 1,1,melri7) = 0. 893 : re( 1,1,nelri7) = 0. 894 : * ip1 895 : re( 2,1,melri7) = alp_d*xr 896 : re( 3,1,melri7) = alp_d*yr 897 : re( 4,1,melri7) = alp_d*zr 898 : re( 2,1,nelri7) = alp_d*xs 899 : re( 3,1,nelri7) = alp_d*ys 900 : re( 4,1,nelri7) = alp_d*zs 901 : * ip2 902 : re( 5,1,melri7) = bet_d*xr 903 : re( 6,1,melri7) = bet_d*yr 904 : re( 7,1,melri7) = bet_d*zr 905 : re( 5,1,nelri7) = bet_d*xs 906 : re( 6,1,nelri7) = bet_d*ys 907 : re( 7,1,nelri7) = bet_d*zs 908 : * ip3 909 : re( 8,1,melri7) = uab_d*xr 910 : re( 9,1,melri7) = uab_d*yr 911 : re(10,1,melri7) = uab_d*zr 912 : re( 8,1,nelri7) = uab_d*xs 913 : re( 9,1,nelri7) = uab_d*ys 914 : re(10,1,nelri7) = uab_d*zs 915 : * ip4 = jp1 916 : re(11,1,melri7) = gam_d*xr 917 : re(12,1,melri7) = gam_d*yr 918 : re(13,1,melri7) = gam_d*zr 919 : re(11,1,nelri7) = gam_d*xs 920 : re(12,1,nelri7) = gam_d*ys 921 : re(13,1,nelri7) = gam_d*zs 922 : * ip5 = jp2 923 : re(14,1,melri7) = del_d*xr 924 : re(15,1,melri7) = del_d*yr 925 : re(16,1,melri7) = del_d*zr 926 : re(14,1,nelri7) = del_d*xs 927 : re(15,1,nelri7) = del_d*ys 928 : re(16,1,nelri7) = del_d*zs 929 : * ip6 = jp3 930 : re(17,1,melri7) = ugd_d*xr 931 : re(18,1,melri7) = ugd_d*yr 932 : re(19,1,melri7) = ugd_d*zr 933 : re(17,1,nelri7) = ugd_d*xs 934 : re(18,1,nelri7) = ugd_d*ys 935 : re(19,1,nelri7) = ugd_d*zs 936 : * on transpose 937 : do 710 ic=2,nligrp 938 : re(1,ic,melri7)=re(ic,1,melri7) 939 : re(1,ic,nelri7)=re(ic,1,nelri7) 940 : 710 continue 941 : * 942 : 700 continue 943 : * 944 : if (nelrig.ne.nelri7) then 945 : nelrig=nelri7 946 : segadj,xmatri 947 : nbelem=nelri7 948 : nbnn =7 949 : nbsous=0 950 : nbref =0 951 : segadj,ipt2 952 : endif 953 : segdes,xmatri,ipt2 954 : endif 955 : 956 : C* SEGDES,mrigid 957 : IPRIG2 = mrigid 958 : 959 : 9910 CONTINUE 960 : SEGSUP,icpr,mwrk 961 : if (mwrkf.ne.0) segsup,mwrkf 962 : 9900 CONTINUE 963 : SEGDES,meleme 964 : 965 : RETURN 966 : END 967 : 968 :
© Cast3M 2003 - All rights reserved.
Disclaimer