1 : C IMPOS2 SOURCE FANDEUR 12/03/15 21:26:27 7314 2 : 3 : * impo bloc en 2D 4 : 5 : SUBROUTINE IMPOS2 6 : 7 : IMPLICIT INTEGER(I-N) 8 : IMPLICIT REAL*8 (A-H,O-Z) 9 : 10 : -INC CCOPTIO 11 : 12 : -INC SMCOORD 13 : -INC SMELEME 14 : -INC SMRIGID 15 : -INC SMCHPOI 16 : 17 : segment mfopa1 18 : integer lpoin(nfopa1,3) 19 : real*8 xpoin(nfopa1,2) 20 : endsegment 21 : 22 : character*4 modepl(4),moforc(4) 23 : 24 : data modepl /'UX ','UY ','UR ','UZ '/ 25 : data moforc /'FX ','FY ','FR ','FZ '/ 26 : * 27 : idimp1 = IDIM + 1 28 : * 29 : * Lecture du maillage support des conditions de contact 30 : * 30 : call lirobj('MAILLAGE',ipt1,1,iretou) 31 : if (ierr.ne.0) return 32 : * 33 : * Activation du maillage et petites verifications 34 : * 35 : segact ipt1 36 : nbno1 = ipt1.num(/1) 37 : nbel1 = ipt1.num(/2) 37 : if (ipt1.lisous(/1).ne.0) call erreur(25) 37 : if (ipt1.itypel.ne.22) call erreur(16) 37 : if (nbno1.ne.4 .and. nbno1.ne.5) call erreur(16) 38 : if (ierr.ne.0) goto 900 39 : * 40 : * Indice des ddls concernes en fonction du mode 2D 41 : * 42 : imo = 1 43 : if (ifour.eq.0) imo = 3 44 : * 45 : * Increment sur le nombre d'elements de contact retenus et utilises 46 : * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les 47 : * longueurs des segments adequats 48 : incnel = 2000 49 : * 50 : * Creation de la raideur des conditions de contact 51 : * Remplissage de l'entete commun 52 : * 53 : nrigel = 1 54 : segini mrigid 55 : ichole = 0 56 : imgeo1 = 0 57 : imgeo2 = 0 58 : isupeq = 0 59 : iforig = ifomod 60 : coerig(1)=1.d0 61 : mtymat='RIGIDITE' 62 : * 63 : * Creation du chpoint de depi 64 : * 65 : nat=1 66 : nsoupo=1 67 : segini mchpoi 68 : mtypoi='DEPIMP' 69 : mochde='engendré par impose' 70 : ifopoi=ifomod 71 : jattri(1)=2 72 : * 73 : nc=1 74 : segini msoupo 75 : ipchp(1)=msoupo 76 : nocomp(1)='FLX ' 77 : nocons(1)=' ' 78 : * 79 : nbnn =1 80 : nbelem=incnel 81 : nbref =0 82 : nbsous=0 83 : segini ipt7 84 : ipt7.itypel=1 85 : igeoc=ipt7 86 : * 87 : n=incnel 88 : segini mpoval 89 : ipoval=mpoval 90 : * 91 : * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval 92 : nelch = 0 93 : * 94 : *======================================================================= 95 : * Formulation "faible" du contact : relation segment-segment (type 0) 96 : *======================================================================= 97 : if (nbno1.ne.5) goto 500 98 : * 99 : * Creation du meleme associe a la relation 100 : * 101 : nbnn =5 102 : nbelem=incnel 103 : nbsous=0 104 : nbref =0 105 : segini meleme 106 : itypel=22 107 : irigel(1,nrigel) = meleme 108 : * 109 : * Creation du descriptif commun a toutes les raideurs 110 : * 111 : nligrp=9 112 : nligrd=nligrp 113 : segini,descr 114 : lisinc(1)='LX ' 115 : lisdua(1)='FLX ' 116 : noelep(1)=1 117 : noeled(1)=1 118 : do 100 i = 2, nligrp, 2 119 : lisinc(i )=modepl(imo) 120 : lisinc(i+1)=modepl(imo+1) 121 : lisdua(i )=moforc(imo) 122 : lisdua(i+1)=moforc(imo+1) 123 : noelep(i )=(i+2)/2 124 : noelep(i+1)=noelep(i) 125 : noeled(i )=noelep(i) 126 : noeled(i+1)=noelep(i) 127 : 100 continue 128 : segdes,descr 129 : irigel(3,nrigel) = descr 130 : * 131 : * creation du segment xmatri 132 : * 133 : nelrig=incnel 134 : segini xmatri 135 : irigel(4,nrigel) = xmatri 136 : * 137 : * ce qu'on cree est unilateral 138 : * 139 : irigel(6,nrigel) = 1 140 : * 141 : * ce qu'on cree est symetrique 142 : * 143 : irigel(7,nrigel) = 0 144 : * 145 : * Nombre d'elements crees dans meleme=irigel(nrigel,1), ipt7 et mpoval 146 : nelri0 = 0 147 : * 148 : * boucle sur le maillage support 149 : * 150 : do 110 iel = 1, nbel1 151 : * 152 : ip1 = ipt1.num(2,iel) 153 : ip2 = ipt1.num(3,iel) 154 : ip3 = ipt1.num(4,iel) 155 : ip4 = ipt1.num(5,iel) 156 : * write(ioimp,*) iel,ip1,ip2,ip3,ip4,ipt1.num(1,iel) 157 : * verification (provisoire) que pas de noeuds doubles 158 : if (ip1.eq.ip3) goto 110 159 : if (ip1.eq.ip4) goto 110 160 : if (ip2.eq.ip3) goto 110 161 : if (ip2.eq.ip4) goto 110 162 : * 163 : ipv = (ip1-1)*idimp1 164 : xp1 = xcoor(ipv+1) 165 : yp1 = xcoor(ipv+2) 166 : ipv = (ip2-1)*idimp1 167 : xp2 = xcoor(ipv+1) 168 : yp2 = xcoor(ipv+2) 169 : xp12 = xp2 - xp1 170 : yp12 = yp2 - yp1 171 : * 172 : ipv = (ip3-1)*idimp1 173 : xp3 = xcoor(ipv+1) 174 : yp3 = xcoor(ipv+2) 175 : ipv = (ip4-1)*idimp1 176 : xp4 = xcoor(ipv+1) 177 : yp4 = xcoor(ipv+2) 178 : xp34 = xp4 - xp3 179 : yp34 = yp4 - yp3 180 : * 181 : * orientations respectives correctes des 2 segments : 182 : * "normales de sens opposes" = produit scalaire negatif ou nul 183 : scal = xp12*xp34 + yp12*yp34 184 : * write(ioimp,*) xp12,yp12,xp34,yp34,scal 185 : if (scal.gt.0.) goto 110 186 : * 187 : * critere d'acceptation de l'élément : 188 : * angles des diagonales 189 : xl13 = sqrt((xp3-xp1)**2+(yp3-yp1)**2) 190 : xl14 = sqrt((xp4-xp1)**2+(yp4-yp1)**2) 191 : xl23 = sqrt((xp3-xp2)**2+(yp3-yp2)**2) 192 : xl24 = sqrt((xp4-xp2)**2+(yp4-yp2)**2) 193 : sca312 = (xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2) 194 : sca412 = (xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2) 195 : sca134 = (xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4) 196 : sca234 = (xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4) 197 : * write(ioimp,*) sca312/(xl13*xl23),sca412/(xl14*xl24), 198 : * & sca134/(xl13*xl14),sca234/(xl23*xl24) 199 : if (sca312/(xl13*xl23).gt.0.7 .and. 200 : > sca412/(xl14*xl24).gt.0.7 .and. 201 : > sca134/(xl13*xl14).gt.0.7 .and. 202 : > sca234/(xl23*xl24).gt.0.7) goto 110 203 : * 204 : * Quelle est la relation ??? 205 : * 206 : * direction du contact unitaire 207 : * 208 : xr = xp12 - xp34 209 : yr = yp12 - yp34 210 : sr = sqrt(xr*xr + yr*yr) 211 : xr = xr/sr 212 : yr = yr/sr 213 : * write(ioimp,*) 'direction contact',xr,yr,yr,-xr 214 : * 215 : * projection des points sur la direction du contact 216 : * 217 : xl1 = xp1*xr + yp1*yr 218 : xl2 = xp2*xr + yp2*yr 219 : xl3 = xp3*xr + yp3*yr 220 : xl4 = xp4*xr + yp4*yr 221 : 222 : * write(ioimp,*) 'projection pts sur contact',xl1,xl2,xl3,xl4 223 : * 224 : * calcul de l'intersection des projections 225 : * 226 : xm1 = min(xl1,xl2) 227 : xm2 = max(xl1,xl2) 228 : xm3 = min(xl3,xl4) 229 : xm4 = max(xl3,xl4) 230 : * write(ioimp,*) ' xmi',xm1,xm2,xm3,xm4 231 : * 232 : * critere de precision sur l'intersection 233 : * 234 : xcr = max(xm2-xm1,xm4-xm3)*1.d-4 235 : * 236 : * taille de l'intersection 237 : * 238 : xi = max(xm1,xm3) 239 : xj = min(xm2,xm4) 240 : xl = xj - xi 241 : * write(ioimp,*) ' intersection',xi,xj,xl,xcr 242 : * 243 : if (xl.le.xcr) goto 110 244 : * 245 : * distance des points a leur projection 246 : * 247 : d1 = (xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr 248 : d2 = (xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr 249 : d3 = (xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr 250 : d4 = (xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr 251 : * 252 : * coordonnées paramétriques de l'intersection sur les segments 1-2 et 3-4 253 : * 254 : xi2 = (xi-xl1) / (xl2-xl1) 255 : xi1 = (xl2-xi) / (xl2-xl1) 256 : 257 : xj2 = (xj-xl1) / (xl2-xl1) 258 : xj1 = (xl2-xj) / (xl2-xl1) 259 : 260 : xi4 = (xi-xl3) / (xl3-xl4) 261 : xi3 = (xl4-xi) / (xl3-xl4) 262 : 263 : xj4 = (xj-xl3) / (xl3-xl4) 264 : xj3 = (xl4-xj) / (xl3-xl4) 265 : * 266 : * write(ioimp,*) ' xi1 xi2 xi3 xi4 xj1 xj2 xj3 xj4 xl ' 267 : * write(ioimp,*) xi1,xi2,xi3,xi4,xj1,xj2,xj3,xj4,xl 268 : * write(ioimp,*) ' d1 d2 d3 d4 ',d1,d2,d3,d4 269 : * 270 : * surface actuelle 271 : * 272 : sc = ((xi1+xj1)*d1+(xi2+xj2)*d2+(xi3+xj3)*d3+(xi4+xj4)*d4)*0.5 273 : * write(ioimp,*) 'Surface actuelle :',sc 274 : * 275 : * on a un element ou imposer la relation a ajouter 276 : * 277 : nelri0 = nelri0 + 1 278 : nelch = nelch +1 279 : * 280 : * on ajuste les differents segments 281 : * 282 : if (nelri0.gt.nelrig) then 283 : nelrig = nelrig + incnel 284 : segadj,xmatri 285 : nbelem = nbelem + incnel 286 : nbnn = 5 287 : segadj,meleme 288 : nbnn = 1 289 : segadj,ipt7 290 : n = n + incnel 291 : segadj,mpoval 292 : endif 293 : * 294 : * on range dans le meleme 295 : * 296 : num(1,nelri0) = ipt1.num(1,iel) 297 : num(2,nelri0) = ip1 298 : num(3,nelri0) = ip2 299 : num(4,nelri0) = ip3 300 : num(5,nelri0) = ip4 301 : * 302 : * on remplit le xmatri 303 : * 304 : * lambda 305 : re(1,1,nelri0)= 0.d0 306 : * ip1 307 : re(2,1,nelri0)= yr * (xi1+xj1) * 0.5 308 : re(3,1,nelri0)= -xr * (xi1+xj1) * 0.5 309 : * ip2 310 : re(4,1,nelri0)= yr * (xi2+xj2) * 0.5 311 : re(5,1,nelri0)= -xr * (xi2+xj2) * 0.5 312 : * ip3 313 : re(6,1,nelri0)= yr * (xi3+xj3) * 0.5 314 : re(7,1,nelri0)= -xr * (xi3+xj3) * 0.5 315 : * ip4 316 : re(8,1,nelri0)= yr * (xi4+xj4) * 0.5 317 : re(9,1,nelri0)= -xr * (xi4+xj4) * 0.5 318 : * 319 : * on transpose 320 : do 120 ic = 2, nligrp 321 : re(1,ic,nelri0)=re(ic,1,nelri0) 322 : 120 continue 323 : * 324 : * Le reste est nul 325 : * 326 : * remplissage du champoint de depi et du maillage support 327 : * 328 : ipt7.num(1,nelch)=ipt1.num(1,iel) 329 : vpocha(nelch,1) = -sc 330 : * 331 : 110 continue 332 : 333 : * write (ioimp,*) ' nb relation type 0 ',nbelem,n,nelri0 334 : 335 : * Ajustement au plus juste puis desactivation des segments lies 336 : * a la rigidite du type 0 337 : if (nelri0.ne.nelrig) then 338 : nelrig = nelri0 339 : segadj,xmatri 340 : nbelem = nelri0 341 : nbnn = 5 342 : segadj,meleme 343 : endif 344 : segdes,meleme,xmatri 345 : * 346 : * S'il n'y a pas d'elements en contact, alors pas de relation unilaterale 347 : * if (nelri0.eq.0) irigel(6,nrigel)=0 348 : * 349 : GOTO 1000 350 : 351 : *======================================================================= 352 : * Formulation "forte" (standard) du contact : 353 : *======================================================================= 354 : 500 continue 355 : * 356 : * Nombre de "cas pathologiques" de type noeud-noeud (type 1) 357 : nelpath = 0 358 : nfopa1 = incnel 359 : segini,mfopa1 360 : * 361 : * Relation type 2 : noeud-segment 362 : *--------------------------------- 363 : * creation du meleme associe a la relation 364 : * 365 : nbnn = 4 366 : nbelem = incnel 367 : nbsous = 0 368 : nbref = 0 369 : segini meleme 370 : itypel=22 371 : irigel(1,nrigel)=meleme 372 : * 373 : * creation du descriptif commun a toutes les raideurs 374 : * 375 : nligrp=7 376 : nligrd=nligrp 377 : segini descr 378 : lisinc(1)='LX ' 379 : lisdua(1)='FLX ' 380 : noelep(1)=1 381 : noeled(1)=1 382 : do 510 i=2,nligrp,2 383 : lisinc(i )=modepl(imo) 384 : lisinc(i+1)=modepl(imo+1) 385 : lisdua(i )=moforc(imo) 386 : lisdua(i+1)=moforc(imo+1) 387 : noelep(i )=(i+2)/2 388 : noelep(i+1)=noelep(i) 389 : noeled(i )=noelep(i) 390 : noeled(i+1)=noelep(i) 391 : 510 continue 392 : segdes,descr 393 : irigel(3,nrigel)=descr 394 : * 395 : * creation du segment xmatri 396 : * 397 : nelrig=incnel 398 : segini xmatri 399 : irigel(4,nrigel)=xmatri 400 : * 401 : * ce qu'on cree est unilateral 402 : * 403 : irigel(6,nrigel)=1 404 : * 405 : * ce qu'on cree est symetrique 406 : * 407 : irigel(7,nrigel)=0 408 : * 409 : * Nombre d'elements dans la rigidite de type 2 410 : nelri2=0 411 : * 412 : * boucle sur le maillage support 413 : * 414 : do 520 iel = 1,nbel1 415 : ip1 = ipt1.num(2,iel) 416 : ip2 = ipt1.num(3,iel) 417 : jp = ipt1.num(4,iel) 418 : *D write(ioimp,*) iel,ip1,ip2,jp 419 : * 420 : * verification que pas relation du point sur lui meme 421 : if (jp.eq.ip1) goto 520 422 : if (jp.eq.ip2) goto 520 423 : 424 : ipv = (ip1-1)*idimp1 425 : xp1 = xcoor(ipv+1) 426 : yp1 = xcoor(ipv+2) 427 : ipv = (ip2-1)*idimp1 428 : xp2 = xcoor(ipv+1) 429 : yp2 = xcoor(ipv+2) 430 : ipv = (jp-1)*idimp1 431 : xp = xcoor(ipv+1) 432 : yp = xcoor(ipv+2) 433 : 434 : x12 = xp2 - xp1 435 : y12 = yp2 - yp1 436 : x1p = xp - xp1 437 : y1p = yp - yp1 438 : x2p = xp - xp2 439 : y2p = yp - yp2 440 : 441 : * distance signee de p a la ligne 1-2 442 : dp12s = x1p * y12 - y1p * x12 443 : * write(ioimp,*) 'dist. signee',dp12s 444 : * verif que le point est du bon cote du segment 445 : * if (dp12s.lt.0.) goto 520 446 : 447 : * verification si autre point dans le cercle de selection 448 : d1p = x1p*x1p + y1p*y1p 449 : d2p = x2p*x2p + y2p*y2p 450 : scal = (x1p*x2p + y1p*y2p) / sqrt(d1p*d2p) 451 : * write(ioimp,*) 'cos(angle_1p2)',scal 452 : if (scal.gt.0.7) goto 520 453 : * 454 : * on a un point ou imposer la relation 455 : * Quelle est la relation ??? 456 : * 457 : * direction de la relation 458 : * 459 : * calcul du pt a mettre en relation avec ip : alpha ip2 + (1-alpha) ip1 460 : * 461 : sr2 = x12*x12 + y12*y12 462 : alpha = (x1p*x12 + y1p*y12) / sr2 463 : * 464 : sr = sqrt(sr2) 465 : xr = x12 / sr 466 : yr = y12 / sr 467 : * traitement des cas pathologiques (on est sur un point du segment) 468 : if (abs(alpha).le.1e-2) then 469 : nelpath = nelpath + 1 470 : if (nelpath.gt.nfopa1) then 471 : nfopa1 = nfopa1 + incnel 472 : segadj,mfopa1 473 : endif 474 : lpoin(nelpath,1)=ip1 475 : lpoin(nelpath,2)=jp 476 : lpoin(nelpath,3)=ipt1.num(1,iel) 477 : xpoin(nelpath,1)= yr 478 : xpoin(nelpath,2)=-xr 479 : * write(ioimp,*) ' cas pathologique 1 (2)',alpha,ip1,jp,yr,-xr 480 : goto 520 481 : else if (abs(1.-alpha).le.1e-2) then 482 : nelpath = nelpath + 1 483 : if (nelpath.gt.nfopa1) then 484 : nfopa1 = nfopa1 + incnel 485 : segadj,mfopa1 486 : endif 487 : lpoin(nelpath,1)=ip2 488 : lpoin(nelpath,2)=jp 489 : lpoin(nelpath,3)=ipt1.num(1,iel) 490 : xpoin(nelpath,1)= yr 491 : xpoin(nelpath,2)=-xr 492 : * write(ioimp,*) ' cas pathologique 1 (1)',alpha,ip1,jp,yr,-xr 493 : goto 520 494 : endif 495 : * 496 : * on a un element ou imposer la relation a ajouter 497 : * 498 : nelri2 = nelri2 + 1 499 : nelch = nelch + 1 500 : * 501 : * on ajuste les differents segments 502 : * 503 : if (nelri2.gt.nelrig) then 504 : nelrig = nelrig + incnel 505 : segadj,xmatri 506 : nbelem = nbelem + incnel 507 : nbnn = 4 508 : segadj,meleme 509 : nbnn = 1 510 : segadj,ipt7 511 : n = n + incnel 512 : segadj,mpoval 513 : endif 514 : * 515 : * on range dans le meleme 516 : * 517 : num(1,nelri2)=ipt1.num(1,iel) 518 : num(2,nelri2)=ip1 519 : num(3,nelri2)=ip2 520 : num(4,nelri2)=jp 521 : * 522 : * on remplit le xmatri 523 : * 524 : * lambda 525 : re(1,1,nelri2)= 0.d0 526 : * ip1 527 : re(2,1,nelri2)= yr * (1.-alpha) 528 : re(3,1,nelri2)= -xr * (1.-alpha) 529 : * ip2 530 : re(4,1,nelri2)= yr * alpha 531 : re(5,1,nelri2)= -xr * alpha 532 : * jp 533 : re(6,1,nelri2)= -yr 534 : re(7,1,nelri2)= xr 535 : * on transpose 536 : re(1,2,nelri2) = re(2,1,nelri2) 537 : re(1,3,nelri2) = re(3,1,nelri2) 538 : re(1,4,nelri2) = re(4,1,nelri2) 539 : re(1,5,nelri2) = re(5,1,nelri2) 540 : re(1,6,nelri2) = re(6,1,nelri2) 541 : re(1,7,nelri2) = re(7,1,nelri2) 542 : * le reste est nul 543 : * 544 : * remplissage du champoint de depi 545 : * 546 : ipt7.num(1,nelch) = ipt1.num(1,iel) 547 : vpocha(nelch,1) = dp12s / sr 548 : *D write(ioimp,*) ' jeu type 2 ',nelri2,nelch,vpocha(nelch,1), 549 : *D & ipt7.num(1,nelch),ip1,ip2,jp 550 : * 551 : 520 continue 552 : * 553 : *D write (ioimp,*) ' nb relation type 2 ',nelri2 554 : * 555 : * Ajustement au plus juste puis desactivation des segments lies 556 : * a la rigidite du type 2 557 : if (nelri2.ne.nelrig) then 558 : nelrig = nelri2 559 : segadj,xmatri 560 : nbelem = nelri2 561 : nbnn = 4 562 : segadj,meleme 563 : endif 564 : segdes,meleme,xmatri 565 : * 566 : * si il n'y a rien on dit que pas unilateral pour pas passer dans unilater 567 : * if (nbelem.eq.0) irigel(6,nrigel) = 0 568 : * 569 : * Relation du type 1 : noeud-noeud 570 : *----------------------------------- 571 : if (nelpath.eq.0) goto 700 572 : * Il y en a au maximum nelpath ! 573 : * 574 : * Ajustement du chpoint et du maillage support associe 575 : * 576 : nbnn = 1 577 : nbelem = nelch + nelpath 578 : nbsous = 0 579 : nbref = 0 580 : segadj,ipt7 581 : n = nelch + nelpath 582 : segadj,mpoval 583 : * 584 : * creation de la raideur de conditions 585 : * 586 : nrigel = irigel(/2)+1 587 : segadj mrigid 588 : coerig(nrigel)=1.d0 589 : * 590 : * creation du meleme associe a la relation 591 : * 592 : nbsous = 0 593 : nbref = 0 594 : nbnn = 3 595 : nbelem = nelpath 596 : segini meleme 597 : itypel=22 598 : irigel(1,nrigel)=meleme 599 : * 600 : * creation des descriptif commun a toutes les raideurs 601 : * 602 : nligrp=5 603 : nligrd=nligrp 604 : segini descr 605 : lisinc(1)='LX ' 606 : lisdua(1)='FLX ' 607 : noelep(1)=1 608 : noeled(1)=1 609 : do 610 i=2, nligrp,2 610 : lisinc(i )=modepl(imo) 611 : lisinc(i+1)=modepl(imo+1) 612 : lisdua(i )=moforc(imo) 613 : lisdua(i+1)=moforc(imo+1) 614 : noelep(i )=(i+2)/2 615 : noelep(i+1)=noelep(i) 616 : noeled(i )=noelep(i) 617 : noeled(i+1)=noelep(i) 618 : 610 continue 619 : segdes,descr 620 : irigel(3,nrigel)=descr 621 : * 622 : * creation du segment imatri 623 : * 624 : nelrig = nelpath 625 : segini xmatri 626 : irigel(4,nrigel)=xmatri 627 : * 628 : * ce qu'on cree est unilateral 629 : * 630 : irigel(6,nrigel)=1 631 : * 632 : * ce qu'on cree est symetrique 633 : * 634 : irigel(7,nrigel)=0 635 : * 636 : nelri1 = 0 637 : * 638 : do 620 iel=1,nelpath 639 : ip1 = lpoin(iel,1) 640 : if (ip1.eq.0) goto 620 641 : *D write (ioimp,*) 'point en relation point point ',ip1 642 : lpoin(iel,1)=0 643 : xp1=xcoor((ip1-1)*idimp1+1) 644 : yp1=xcoor((ip1-1)*idimp1+2) 645 : jp =lpoin(iel,2) 646 : xp =xcoor((jp -1)*idimp1+1) 647 : yp =xcoor((jp -1)*idimp1+2) 648 : il1=lpoin(iel,3) 649 : xn = xpoin(iel,1) 650 : yn = xpoin(iel,2) 651 : sn1 = sqrt(xn*xn+yn*yn) 652 : * write (6,*) ' xn yn avant 630 ',xn,yn 653 : do 630 jel=iel+1,nelpath 654 : if (lpoin(jel,1).ne.ip1.or.lpoin(jel,2).ne.jp) goto 630 655 : * write (ioimp,*) 'point en double ' 656 : lpoin(jel,1)=0 657 : xn=xn+xpoin(jel,1) 658 : yn=yn+xpoin(jel,2) 659 : 630 continue 660 : do 640 jel=iel+1,nelpath 661 : if (lpoin(jel,2).ne.ip1.or.lpoin(jel,1).ne.jp) goto 640 662 : * write (6,*) 'point en double ' 663 : lpoin(jel,1)=0 664 : xn=xn-xpoin(jel,1) 665 : yn=yn-xpoin(jel,2) 666 : 640 continue 667 : *D write(ioimp,*) ' xn yn sn ',xn,yn,sn 668 : sn=sqrt(xn*xn+yn*yn) 669 : if (sn.lt.sn1*1.e-3) then 669 : call erreur(942) 670 : return 671 : endif 672 : xn=xn/sn 673 : yn=yn/sn 674 : *D write (ioimp,*) ' xn yn sn apres 640 ',xn,yn,sn,ip1,jp 675 : * 676 : nelri1 = nelri1 + 1 677 : nelch = nelch + 1 678 : * 679 : * on range dans le meleme 680 : * 681 : num(1,nelri1)=il1 682 : num(2,nelri1)=ip1 683 : num(3,nelri1)=jp 684 : * 685 : * on remplit le xmatri 686 : * 687 : * lambda 688 : re(1,1,nelri1)= 0. 689 : * ip1 690 : re(2,1,nelri1)= xn 691 : re(3,1,nelri1)= yn 692 : * jp 693 : re(4,1,nelri1)= -xn 694 : re(5,1,nelri1)= -yn 695 : * on transpose 696 : re(1,2,nelri1) = re(2,1,nelri1) 697 : re(1,3,nelri1) = re(3,1,nelri1) 698 : re(1,4,nelri1) = re(4,1,nelri1) 699 : re(1,5,nelri1) = re(5,1,nelri1) 700 : * le reste est nul 701 : * 702 : * remplissage du champoint de depi 703 : * 704 : ipt7.num(1,nelch) = il1 705 : scal=(xp-xp1)*xn+(yp-yp1)*yn 706 : vpocha(nelch,1) = scal 707 : * write(ioimp,*) ' jeu type 1 ',scal,ip1,jp,il1,il2 708 : * 709 : 620 continue 710 : * 711 : * write (ioimp,*) ' nb relation type 1 ',nelri1 712 : * 713 : * Ajustement au plus juste puis desactivation des segments lies 714 : * a la rigidite du type 1 715 : if (nelri1.ne.nelrig) then 716 : nelrig = nelri1 717 : segadj,xmatri 718 : nbelem = nelri1 719 : nbnn = 3 720 : segadj,meleme 721 : endif 722 : segdes,meleme,xmatri 723 : * 724 : * if (nelri1.eq.0) irigel(6,nrigel)=0 725 : 726 : 700 CONTINUE 727 : * 728 : segsup,mfopa1 729 : * GOTO 1000 730 : * 731 : *---------------------------- 732 : * on renvoie le resultat 733 : *---------------------------- 734 : 1000 CONTINUE 735 : * 736 : * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7 737 : * puis desactivation du chpoint 738 : if (n.ne.nelch) then 739 : n = nelch 740 : segadj,mpoval 741 : nbnn = 1 742 : nbelem = nelch 743 : nbsous = 0 744 : nbref = 0 745 : segadj,ipt7 746 : endif 747 : segdes,mpoval,ipt7 748 : segdes,msoupo,mchpoi 749 : * Desctivation de la matrice de raideur de contact 750 : segdes,mrigid 751 : * 752 : * Ecriture des objets resultats 752 : call ecrobj('RIGIDITE',mrigid) 752 : call ecrobj('CHPOINT',mchpoi) 753 : 754 : 900 continue 755 : segdes,ipt1 756 : 757 : return 758 : end 759 : 760 : 761 : 762 :
© Cast3M 2003 - All rights reserved.
Disclaimer