1 : C IMPO32 SOURCE PV 13/12/13 21:15:04 7883 2 : 3 : * impo bloc en 3D 4 : 5 : SUBROUTINE IMPO32 6 : 7 : IMPLICIT INTEGER(I-N) 8 : IMPLICIT REAL*8 (A-H,O-Z) 9 : 10 : -INC CCOPTIO 11 : -INC CCREEL 12 : -INC CCGEOME 13 : 14 : -INC SMCOORD 15 : -INC SMELEME 16 : -INC SMRIGID 17 : -INC SMCHPOI 18 : 19 : * Contact en formulation forte : 20 : * Segments utiles dans les cas dits pathologiques 21 : segment mfopa1 22 : integer lpoin(nfopa1,3) 23 : real*8 xpoin(nfopa1,3) 24 : endsegment 25 : segment mfopa2 26 : integer lsegm(nfopa2,4) 27 : real*8 xsegm(nfopa2,3) 28 : endsegment 29 : * directions moyennes aux sommets 30 : segment mfopa3 31 : real*8 xms(nbpts),yms(nbpts),zms(nbpts) 32 : endsegment 33 : 34 : * Contact en formulation faible : 35 : * Segment utile pour le calcul de l'intersection des 2 triangles 36 : * projetes dans le plan de contact intermediaire 37 : segment mfaible 38 : real*8 cPT1C(3,3),cPT2C(3,3), cPT1A(4,2),cPT2A(4,2) 39 : real*8 cPIn0(6,2), cPIn(6,2), test(6) 40 : real*8 vT1A(3,2), vT2A(3,2) 41 : real*8 SuT1A, SuT2A, SuPIn 42 : *DBG-F real*8 xGIn,yGIn, b1T1,b2T1,b3T1, b1T2,b2T2,b3T2 43 : endsegment 44 : 45 : PARAMETER ( X1s3=0.333333333333333333333333333333333333333333D0 ) 46 : 47 : character*4 modepl(3),moforc(3) 48 : 49 : data modepl /'UX ','UY ','UZ '/ 50 : data moforc /'FX ','FY ','FZ '/ 51 : * 52 : idimp1 = IDIM + 1 53 : * 54 : * Lecture du maillage support des conditions de contact 55 : * 55 : call lirobj('MAILLAGE',ipt1,1,iretou) 56 : if (ierr.ne.0) return 57 : * 58 : * Activation du maillage et petites verifications 59 : * 60 : segact ipt1 61 : nbno1 = ipt1.num(/1) 62 : nbel1 = ipt1.num(/2) 62 : if (ipt1.lisous(/1).ne.0) call erreur(25) 62 : if (ipt1.itypel.ne.22) call erreur(16) 62 : if (nbno1.ne.5 .and. nbno1.ne.7) call erreur(16) 63 : if (ierr.ne.0) goto 900 64 : * 65 : * Increment sur le nombre d'elements de contact retenus et utilises 66 : * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les 67 : * longueurs des segments adequats 68 : incnel = 2000 69 : * 70 : * Creation de la raideur des conditions de contact 71 : * Remplissage de l'entete commun 72 : * 73 : nrigel = 1 74 : segini,mrigid 75 : ichole = 0 76 : imgeo1 = 0 77 : imgeo2 = 0 78 : isupeq = 0 79 : iforig = ifomod 80 : coerig(1) = 1. 81 : mtymat='RIGIDITE' 82 : * 83 : * Creation du chpoint de depi 84 : * 85 : nat=1 86 : nsoupo=1 87 : segini mchpoi 88 : mtypoi='DEPIMP' 89 : mochde='engendré par impose' 90 : ifopoi=ifomod 91 : jattri(1)=2 92 : * 93 : nc=1 94 : segini msoupo 95 : ipchp(1)=msoupo 96 : nocomp(1)='FLX ' 97 : nocons(1)=' ' 98 : * 99 : nbnn =1 100 : nbelem=incnel 101 : nbref =0 102 : nbsous=0 103 : segini ipt7 104 : ipt7.itypel=1 105 : igeoc=ipt7 106 : * 107 : n=incnel 108 : segini mpoval 109 : ipoval=mpoval 110 : * 111 : * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval 112 : nelch = 0 113 : 114 : *======================================================================= 115 : * Formulation "forte" (standard) du contact : 116 : *======================================================================= 117 : * Element de contact 3D a 5 noeuds 118 : if (nbno1.eq.7) goto 500 119 : * Calcul des directions moyennes 120 : nbpts = xcoor(/1)/(idim+1) 121 : segini mfopa3 122 : do 820 iel=1,nbel1 123 : * 124 : ip1=ipt1.num(2,iel) 125 : ip2=ipt1.num(3,iel) 126 : ip3=ipt1.num(4,iel) 127 : ipv1 = (ip1-1)*idimp1 128 : xp1 = xcoor(ipv1+1) 129 : yp1 = xcoor(ipv1+2) 130 : zp1 = xcoor(ipv1+3) 131 : ipv2 = (ip2-1)*idimp1 132 : xp2 = xcoor(ipv2+1) 133 : yp2 = xcoor(ipv2+2) 134 : zp2 = xcoor(ipv2+3) 135 : ipv3 = (ip3-1)*idimp1 136 : xp3 = xcoor(ipv3+1) 137 : yp3 = xcoor(ipv3+2) 138 : zp3 = xcoor(ipv3+3) 139 : * 140 : * normale au plan (123) 141 : * 142 : x12 = xp2 - xp1 143 : y12 = yp2 - yp1 144 : z12 = zp2 - zp1 145 : x23 = xp3 - xp2 146 : y23 = yp3 - yp2 147 : z23 = zp3 - zp2 148 : xn = y12*z23 - z12*y23 149 : yn = z12*x23 - x12*z23 150 : zn = x12*y23 - y12*x23 151 : sn = sqrt (xn*xn + yn*yn + zn*zn) 152 : xn = xn/sn 153 : yn = yn/sn 154 : zn = zn/sn 155 : xms(ip1)=xms(ip1)+xn 156 : yms(ip1)=yms(ip1)+yn 157 : zms(ip1)=zms(ip1)+zn 158 : xms(ip2)=xms(ip2)+xn 159 : yms(ip2)=yms(ip2)+yn 160 : zms(ip2)=zms(ip2)+zn 161 : xms(ip3)=xms(ip3)+xn 162 : yms(ip3)=yms(ip3)+yn 163 : zms(ip3)=zms(ip3)+zn 164 : 820 continue 165 : do 821 ip = 1,nbpts 166 : sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)+zms(ip)*zms(ip) 167 : sn = max(sqrt(abs(sn)),1d-35) 168 : xms(ip)=xms(ip)/sn 169 : yms(ip)=yms(ip)/sn 170 : zms(ip)=zms(ip)/sn 171 : 821 continue 172 : * 173 : * Nombre de "cas pathologiques" de type noeud-noeud (type 1) 174 : nelpa1 = 0 175 : * Nombre de "cas pathologiques" de type noeud-segment (type 2) 176 : nelpa2 = 0 177 : * Critere de detection des cas pathologiques 178 : dcrit = 1.d-2 179 : * Segments lies aux cas pathologiques (types 1 et 2) 180 : nfopa1 = incnel 181 : nfopa2 = incnel 182 : segini,mfopa1,mfopa2 183 : * 184 : * Relation du type 3 : noeud-triangle 185 : *------------------------------------- 186 : * 187 : * creation du meleme associe a la relation 188 : * 189 : nbsous=0 190 : nbref =0 191 : nbnn =5 192 : nbelem=incnel 193 : segini meleme 194 : itypel=22 195 : irigel(1,nrigel)=meleme 196 : * 197 : * creation du descriptif commun a toutes les raideurs 198 : * 199 : nligrp = 13 200 : nligrd = nligrp 201 : segini,descr 202 : lisinc(1)='LX ' 203 : lisdua(1)='FLX ' 204 : noelep(1)=1 205 : noeled(1)=1 206 : do 10 i=2, nligrp, 3 207 : lisinc(i) =modepl(1) 208 : lisinc(i+1)=modepl(2) 209 : lisinc(i+2)=modepl(3) 210 : lisdua(i) =moforc(1) 211 : lisdua(i+1)=moforc(2) 212 : lisdua(i+2)=moforc(3) 213 : noelep(i) =(i+5)/3 214 : noelep(i+1)=noelep(i) 215 : noelep(i+2)=noelep(i) 216 : noeled(i) =noelep(i) 217 : noeled(i+1)=noelep(i) 218 : noeled(i+2)=noelep(i) 219 : 10 continue 220 : segdes,descr 221 : irigel(3,nrigel) = descr 222 : * 223 : * creation du segment xmatri 224 : * 225 : nelrig = incnel 226 : segini,xmatri 227 : irigel(4,nrigel) = xmatri 228 : * 229 : * ce qu'on cree est unilateral 230 : * 231 : irigel(6,nrigel)=1 232 : * 233 : * ce qu'on cree est symetrique 234 : * 235 : irigel(7,nrigel)=0 236 : * 237 : * Nombre d'elements dans la rigidite du type 3 238 : nelri3 = 0 239 : 240 : * Boucle sur le maillage support du contact(frottement) 241 : * 242 : ip1o = 0 243 : ip2o = 0 244 : ip3o = 0 245 : DO 20 iel=1,nbel1 246 : * 247 : ip1=ipt1.num(2,iel) 248 : ip2=ipt1.num(3,iel) 249 : ip3=ipt1.num(4,iel) 250 : if (ip1.ne.ip1o.or.ip2.ne.ip2o.or.ip3.ne.ip3o) then 251 : ip1o=ip1 252 : ip2o=ip2 253 : ip3o=ip3 254 : * Recuperation des coordonees des noeuds du triangle 255 : ipv = (ip1-1)*idimp1 256 : xp1 = xcoor(ipv+1) 257 : yp1 = xcoor(ipv+2) 258 : zp1 = xcoor(ipv+3) 259 : ipv = (ip2-1)*idimp1 260 : xp2 = xcoor(ipv+1) 261 : yp2 = xcoor(ipv+2) 262 : zp2 = xcoor(ipv+3) 263 : ipv = (ip3-1)*idimp1 264 : xp3 = xcoor(ipv+1) 265 : yp3 = xcoor(ipv+2) 266 : zp3 = xcoor(ipv+3) 267 : * Centre de gravite du triangle 268 : xg = (xp1+xp2+xp3) * X1s3 269 : yg = (yp1+yp2+yp3) * X1s3 270 : zg = (zp1+zp2+zp3) * X1s3 271 : endif 272 : jp =ipt1.num(5,iel) 273 : ipv = (jp-1)*idimp1 274 : xp = xcoor(ipv+1) 275 : yp = xcoor(ipv+2) 276 : zp = xcoor(ipv+3) 277 : * Verification que pas relation du point sur lui meme 278 : if (jp.eq.ip1) goto 20 279 : if (jp.eq.ip2) goto 20 280 : if (jp.eq.ip3) goto 20 281 : * 282 : * Verification si autre point dans la zone de selection 283 : * verif distance au centre de gravite 284 : d1 = (xg-xp1)**2 + (yg-yp1)**2 + (zg-zp1)**2 285 : d2 = (xg-xp2)**2 + (yg-yp2)**2 + (zg-zp2)**2 286 : d3 = (xg-xp3)**2 + (yg-yp3)**2 + (zg-zp3)**2 287 : dist = max(d1,d2,d3) 288 : ** dp = ((xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2) * 0.25 289 : dp = ((xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2) * 0.64 290 : c*OF A voir -> dp = (xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2 291 : if (dp.gt.dist*2) goto 20 292 : C*DBG write(ioimp,*) 'contact test distance ok',dp,d1,d2,d3 293 : dist = SQRT(dist) 294 : * verif position par rapport aux cotes 295 : * cote 1 2 296 : x12 = xp2 - xp1 297 : y12 = yp2 - yp1 298 : z12 = zp2 - zp1 299 : dv = sqrt((x12*x12)+(y12*y12)+(z12*z12)) 300 : if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10)) 301 : & write(ioimp,*) 'Prob 1.1 - impo32' 302 : if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 1.2 - impo32' 303 : xv = x12 / dv 304 : yv = y12 / dv 305 : zv = z12 / dv 306 : scal = (xp-xp1)*xv + (yp-yp1)*yv + (zp-zp1)*zv 307 : xm = xp1 + scal*xv 308 : ym = yp1 + scal*yv 309 : zm = zp1 + scal*zv 310 : scal = (xg-xm)*xv + (yg-ym)*yv + (zg-zm)*zv 311 : xn = (xg-xm) - scal*xv 312 : yn = (yg-ym) - scal*yv 313 : zn = (zg-zm) - scal*zv 314 : dn = sqrt(xn*xn+yn*yn+zn*zn) 315 : dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) 316 : scal = (xp-xm)*xn + (yp-ym)*yn + (zp-zm)*zn 317 : if (dpm.gt.1d-3*dv .and. scal.lt.-0.8d0*dn*dpm) goto 20 318 : C*DBG write(ioimp,*) 'contact test position 1 ok', 319 : C*DBG & scal/(dn*dpm),scal,dn,dpm 320 : * cote 2 3 321 : x23 = xp3 - xp2 322 : y23 = yp3 - yp2 323 : z23 = zp3 - zp2 324 : dv = sqrt((x23*x23)+(y23*y23)+(z23*z23)) 325 : if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10)) 326 : & write(ioimp,*) 'Prob 2.1 - impo32' 327 : if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 2.2 - impo32' 328 : xv = x23 / dv 329 : yv = y23 / dv 330 : zv = z23 / dv 331 : scal = (xp-xp2)*xv + (yp-yp2)*yv + (zp-zp2)*zv 332 : xm = xp2 + scal*xv 333 : ym = yp2 + scal*yv 334 : zm = zp2 + scal*zv 335 : scal = (xg-xm)*xv + (yg-ym)*yv + (zg-zm)*zv 336 : xn = (xg-xm) - scal*xv 337 : yn = (yg-ym) - scal*yv 338 : zn = (zg-zm) - scal*zv 339 : dn = sqrt(xn*xn+yn*yn+zn*zn) 340 : dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) 341 : scal = (xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn 342 : if (dpm.gt.1d-3*dv .and. scal.lt.-0.8d0*dn*dpm) goto 20 343 : C*DBG write(ioimp,*) 'contact test position 2 ok', 344 : C*DBG & scal/(dn*dpm),scal,dn,dpm 345 : * cote 3 1 346 : x31 = xp1 - xp3 347 : y31 = yp1 - yp3 348 : z31 = zp1 - zp3 349 : dv = sqrt((x31*x31)+(y31*y31)+(z31*z31)) 350 : if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10)) 351 : & write(ioimp,*) 'Prob 3.1 - impo32' 352 : if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 3.2 - impo32' 353 : xv = x31 / dv 354 : yv = y31 / dv 355 : zv = z31 / dv 356 : scal = (xp-xp3)*xv + (yp-yp3)*yv + (zp-zp3)*zv 357 : xm = xp3 + scal*xv 358 : ym = yp3 + scal*yv 359 : zm = zp3 + scal*zv 360 : scal = (xg-xm)*xv + (yg-ym)*yv + (zg-zm)*zv 361 : xn = (xg-xm) - scal*xv 362 : yn = (yg-ym) - scal*yv 363 : zn = (zg-zm) - scal*zv 364 : dn = sqrt(xn*xn+yn*yn+zn*zn) 365 : dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) 366 : scal = (xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn 367 : if (dpm.gt.1d-3*dv .and .scal.lt.-0.8d0*dn*dpm) goto 20 368 : C*DBG write(ioimp,*) 'contact test position 3 ok', 369 : C*DBG & scal/(dn*dpm),scal,dn,dpm 370 : * 371 : * on a un point ou imposer la relation 372 : C*DBG write(ioimp,*) 'contact potentiel ' 373 : * Quelle est la relation ??? 374 : * 375 : * direction de la relation = normale au plan (123) 376 : * 377 : xn = y12*z23 - z12*y23 378 : yn = z12*x23 - x12*z23 379 : zn = x12*y23 - y12*x23 380 : dn = sqrt(xn*xn+yn*yn+zn*zn) 381 : if (.not.(dn.lt.1d-10).and. .not.(dn.gt.1d-10)) 382 : & write(ioimp,*) 'Prob 4.1 - impo32' 383 : if (abs(dn).le.xpetit) write(ioimp,*) 'Prob 4.2 - impo32' 384 : xn = xn / dn 385 : yn = yn / dn 386 : zn = zn / dn 387 : if (abs(xn).lt.1d-12) xn = 0.d0 388 : if (abs(yn).lt.1d-12) yn = 0.d0 389 : if (abs(zn).lt.1d-12) zn = 0.d0 390 : * 391 : * Distance (jeu) du point jp au plan de contact 392 : * Puis calcul de la projection sur le plan de contact 393 : * 394 : xjeu = (xp-xg)*xn + (yp-yg)*yn + (zp-zg)*zn 395 : * verif bon cote (a la tolerance pres) 396 : * write (6,*) ' xjeu ',xjeu,dist 397 : 398 : if (xjeu.lt.-0.8d0*dist) goto 20 399 : 400 : * 401 : * Recherche de ses coordonnées barycentriques 402 : * qui sont les surfaces des sous triangles 403 : * 404 : xb1 = (yp-yp2)*(zp-zp3)-(zp-zp2)*(yp-yp3) 405 : yb1 = (zp-zp2)*(xp-xp3)-(xp-xp2)*(zp-zp3) 406 : zb1 = (xp-xp2)*(yp-yp3)-(yp-yp2)*(xp-xp3) 407 : xb2 = (yp-yp3)*(zp-zp1)-(zp-zp3)*(yp-yp1) 408 : yb2 = (zp-zp3)*(xp-xp1)-(xp-xp3)*(zp-zp1) 409 : zb2 = (xp-xp3)*(yp-yp1)-(yp-yp3)*(xp-xp1) 410 : xb3 = (yp-yp1)*(zp-zp2)-(zp-zp1)*(yp-yp2) 411 : yb3 = (zp-zp1)*(xp-xp2)-(xp-xp1)*(zp-zp2) 412 : zb3 = (xp-xp1)*(yp-yp2)-(yp-yp1)*(xp-xp2) 413 : b1 = xb1*xn + yb1*yn + zb1*zn 414 : b2 = xb2*xn + yb2*yn + zb2*zn 415 : b3 = xb3*xn + yb3*yn + zb3*zn 416 : bt = b1 + b2 + b3 417 : if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10)) 418 : & write(ioimp,*) 'Prob 5.1 - impo32' 419 : if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 5.2 - impo32' 420 : b1 = b1 / bt 421 : b2 = b2 / bt 422 : b3 = b3 / bt 423 : b1c = b1 424 : b2c = b2 425 : b3c = b3 426 : * recalcul de la direction et du jeu en fonction des normales aux sommets 427 : if (b3c.lt.0) then 428 : b3c = 0.d0 429 : b1c = (xp-xp2)*(xp1-xp2)+(yp-yp2)*(yp1-yp2)+(zp-zp2)*(zp1-zp2) 430 : b2c = (xp-xp1)*(xp2-xp1)+(yp-yp1)*(yp2-yp1)+(zp-zp1)*(zp2-zp1) 431 : elseif (b1c.lt.0.d0) then 432 : b1c = 0.d0 433 : b2c = (xp-xp3)*(xp2-xp3)+(yp-yp3)*(yp2-yp3)+(zp-zp3)*(zp2-zp3) 434 : b3c = (xp-xp2)*(xp3-xp2)+(yp-yp2)*(yp3-yp2)+(zp-zp2)*(zp3-zp2) 435 : elseif (b2c.lt.0.d0) then 436 : b2c=0.d0 437 : b3c = (xp-xp1)*(xp3-xp1)+(yp-yp1)*(yp3-yp1)+(zp-zp1)*(zp3-zp1) 438 : b1c = (xp-xp3)*(xp1-xp3)+(yp-yp3)*(yp1-yp3)+(zp-zp3)*(zp1-zp3) 439 : endif 440 : b1c = max(b1c,0.d0) 441 : b2c = max(b2c,0.d0) 442 : b3c = max(b3c,0.d0) 443 : btc = b1c + b2c + b3c 444 : b1c = b1c / btc 445 : b2c = b2c / btc 446 : b3c = b3c / btc 447 : xn = b1c*xms(ip1)+b2c*xms(ip2)+b3c*xms(ip3) 448 : yn = b1c*yms(ip1)+b2c*yms(ip2)+b3c*yms(ip3) 449 : zn = b1c*zms(ip1)+b2c*zms(ip2)+b3c*zms(ip3) 450 : sn = sqrt (xn*xn + yn*yn + zn*zn) 451 : xn=xn/sn 452 : yn=yn/sn 453 : zn=zn/sn 454 : * verif compatibilite avec la normale au poin impactant 455 : if (xms(jp)*xn+yms(jp)*yn+zms(jp)*zn.gt.0.d0) goto 20 456 : xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn 457 : xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn 458 : xjeu3 = (xp-xp3)*xn + (yp-yp3)*yn + (zp-zp3)*zn 459 : xjeu = xjeu1 * b1c + xjeu2 * b2c + xjeu3 * b3c 460 : ** write (6,*) ' normale ',xn,yn,zn,' jeu ',xjeu1,xjeu2,xjeu3 461 : 462 : C*DBG write(ioimp,*) ' b1 b2 b3 ',b1,b2,b3 463 : * creation du support du multiplicateur de lagrange si necessaire 464 : nbpts = xcoor(/1)/(idim+1) 465 : if (ipt1.num(1,iel).eq.ilgni) then 466 : nbpts = nbpts + 1 467 : segadj mcoord 468 : segact ipt1*mod 469 : ipt1.num(1,iel)=nbpts 470 : endif 471 : segact mcoord*mod 472 : xcoor(nbpts*(idim+1)) = dist 473 : * Detection et traitement cas pathologiques 474 : if (b1.lt.dcrit) then 475 : if (b2.lt.dcrit) then 476 : C*DBG write(ioimp,*) 'path 1 (12)',ip3,jp 477 : nelpa1 = nelpa1 + 1 478 : if (nelpa1.gt.nfopa1) then 479 : nfopa1 = nfopa1 + incnel 480 : segadj,mfopa1 481 : endif 482 : lpoin(nelpa1,1)=ip3 483 : lpoin(nelpa1,2)=jp 484 : lpoin(nelpa1,3)=ipt1.num(1,iel) 485 : xpoin(nelpa1,1)=xn 486 : xpoin(nelpa1,2)=yn 487 : xpoin(nelpa1,3)=zn 488 : else if (b3.lt.dcrit) then 489 : C*DBG write (ioimp,*) 'path 1 (13)',ip2,jp 490 : nelpa1 = nelpa1 + 1 491 : if (nelpa1.gt.nfopa1) then 492 : nfopa1 = nfopa1 + incnel 493 : segadj,mfopa1 494 : endif 495 : lpoin(nelpa1,1)=ip2 496 : lpoin(nelpa1,2)=jp 497 : lpoin(nelpa1,3)=ipt1.num(1,iel) 498 : xpoin(nelpa1,1)=xn 499 : xpoin(nelpa1,2)=yn 500 : xpoin(nelpa1,3)=zn 501 : else 502 : C*DBG write (ioimp,*) 'path 2 (1)',ip2,ip3,jp 503 : nelpa2 = nelpa2 + 1 504 : if (nelpa2.gt.nfopa2) then 505 : nfopa2 = nfopa2 + incnel 506 : segadj,mfopa2 507 : endif 508 : lsegm(nelpa2,1)=min(ip2,ip3) 509 : lsegm(nelpa2,2)=max(ip2,ip3) 510 : lsegm(nelpa2,3)=jp 511 : lsegm(nelpa2,4)=ipt1.num(1,iel) 512 : xsegm(nelpa2,1)=xn 513 : xsegm(nelpa2,2)=yn 514 : xsegm(nelpa2,3)=zn 515 : endif 516 : goto 20 517 : else if (b2.lt.dcrit) then 518 : if (b3.lt.dcrit) then 519 : C*DBG write(ioimp,*) 'path 1 (12)',ip3,jp 520 : nelpa1 = nelpa1 + 1 521 : if (nelpa1.gt.nfopa1) then 522 : nfopa1 = nfopa1 + incnel 523 : segadj,mfopa1 524 : endif 525 : lpoin(nelpa1,1)=ip1 526 : lpoin(nelpa1,2)=jp 527 : lpoin(nelpa1,3)=ipt1.num(1,iel) 528 : xpoin(nelpa1,1)=xn 529 : xpoin(nelpa1,2)=yn 530 : xpoin(nelpa1,3)=zn 531 : else 532 : C*DBG write (ioimp,*) 'path 2 (2)',ip1,ip3,jp 533 : nelpa2 = nelpa2 + 1 534 : if (nelpa2.gt.nfopa2) then 535 : nfopa2 = nfopa2 + incnel 536 : segadj,mfopa2 537 : endif 538 : lsegm(nelpa2,1)=min(ip1,ip3) 539 : lsegm(nelpa2,2)=max(ip1,ip3) 540 : lsegm(nelpa2,3)=jp 541 : lsegm(nelpa2,4)=ipt1.num(1,iel) 542 : xsegm(nelpa2,1)=xn 543 : xsegm(nelpa2,2)=yn 544 : xsegm(nelpa2,3)=zn 545 : endif 546 : goto 20 547 : else if (b3.lt.dcrit) then 548 : C*DBG write (ioimp,*) 'path 2 (3)',ip2,ip3,jp 549 : nelpa2 = nelpa2 + 1 550 : if (nelpa2.gt.nfopa2) then 551 : nfopa2 = nfopa2 + incnel 552 : segadj,mfopa2 553 : endif 554 : lsegm(nelpa2,1)=min(ip1,ip2) 555 : lsegm(nelpa2,2)=max(ip1,ip2) 556 : lsegm(nelpa2,3)=jp 557 : lsegm(nelpa2,4)=ipt1.num(1,iel) 558 : xsegm(nelpa2,1)=xn 559 : xsegm(nelpa2,2)=yn 560 : xsegm(nelpa2,3)=zn 561 : goto 20 562 : endif 563 : * 564 : * Ajout d'une relation noeud-triangle 565 : * 566 : nelri3 = nelri3 + 1 567 : nelch = nelch + 1 568 : * 569 : * on ajuste les differents segments si necesaire 570 : * 571 : if (nelri3.gt.nelrig) then 572 : nelrig = nelrig + incnel 573 : segadj,xmatri 574 : nbelem = nbelem + incnel 575 : nbnn = 5 576 : segadj,meleme 577 : nbnn = 1 578 : segadj,ipt7 579 : n = n + incnel 580 : segadj,mpoval 581 : endif 582 : * 583 : * Mise a jour du meleme 584 : num(1,nelri3) = ipt1.num(1,iel) 585 : num(2,nelri3) = ip1 586 : num(3,nelri3) = ip2 587 : num(4,nelri3) = ip3 588 : num(5,nelri3) = jp 589 : * 590 : * Mise a jour de xmatri 591 : * lambda 592 : re( 1,1,nelri3) = 0d0 593 : * ip1 594 : re( 2,1,nelri3) = xn * b1 595 : re( 3,1,nelri3) = yn * b1 596 : re( 4,1,nelri3) = zn * b1 597 : * ip2 598 : re( 5,1,nelri3) = xn * b2 599 : re( 6,1,nelri3) = yn * b2 600 : re( 7,1,nelri3) = zn * b2 601 : * ip3 602 : re( 8,1,nelri3) = xn * b3 603 : re( 9,1,nelri3) = yn * b3 604 : re(10,1,nelri3) = zn * b3 605 : * jp 606 : re(11,1,nelri3) = -xn 607 : re(12,1,nelri3) = -yn 608 : re(13,1,nelri3) = -zn 609 : * on transpose 610 : do 30 ic = 2, nligrp 611 : re(1,ic,nelri3) = re(ic,1,nelri3) 612 : 30 continue 613 : * le reste est nul 614 : * 615 : * remplissage du champoint de depi (jeu) 616 : * 617 : ipt7.num(1,nelch) = ipt1.num(1,iel) 618 : vpocha(nelch,1) = xjeu 619 : C*DBG write(ioimp,*) ' jeu ',xjeu,ip1,ip2,p3,jp,ipt1.num(1,iel) 620 : * 621 : 20 CONTINUE 622 : * 623 : * Ajustement au plus juste puis desactivation des segments lies 624 : * a la rigidite du type 3 625 : if (nelri3.ne.nelrig) then 626 : nelrig = nelri3 627 : segadj,xmatri 628 : nbelem = nelri3 629 : nbnn = 5 630 : segadj,meleme 631 : endif 632 : segdes,meleme,xmatri 633 : * 634 : write(ioimp,*) ' nb relation type 3 ',nelri3 635 : C*DBG if (nelri3.eq.0) irigel(6,nrigel)=0 636 : * 637 : * Relation du type 2 : noeud-segment 638 : *----------------------------------- 639 : if (nelpa2.eq.0) goto 200 640 : * 641 : * Il y a au maximum nelpa2 relations de type 2 ! 642 : * 643 : * Ajustement du chpoint de depi (jeu) : mpoval et ipt7 644 : * 645 : n = nelch + nelpa2 646 : segadj,mpoval 647 : nbnn = 1 648 : nbelem = nelch + nelpa2 649 : nbsous = 0 650 : nbref = 0 651 : segadj,ipt7 652 : * 653 : * creation de la raideur de conditions de type 2 654 : * 655 : nrigel = nrigel+1 656 : segadj,mrigid 657 : coerig(nrigel)=1.d0 658 : * 659 : * creation du meleme associe a la relation 660 : * 661 : nbnn = 4 662 : nbelem = nelpa2 663 : nbsous = 0 664 : nbref = 0 665 : segini,meleme 666 : itypel=22 667 : irigel(1,nrigel)=meleme 668 : * 669 : * creation du descriptif commun a toutes les raideurs 670 : * 671 : nligrp=10 672 : nligrd=nligrp 673 : segini descr 674 : lisinc(1)='LX' 675 : lisdua(1)='FLX' 676 : noelep(1)=1 677 : noeled(1)=1 678 : do 110 i = 2, 10, 3 679 : lisinc(i )=modepl(1) 680 : lisinc(i+1)=modepl(2) 681 : lisinc(i+2)=modepl(3) 682 : lisdua(i )=moforc(1) 683 : lisdua(i+1)=moforc(2) 684 : lisdua(i+2)=moforc(3) 685 : noelep(i )=(i+5)/3 686 : noelep(i+1)=noelep(i) 687 : noelep(i+2)=noelep(i) 688 : noeled(i )=noelep(i) 689 : noeled(i+1)=noelep(i) 690 : noeled(i+2)=noelep(i) 691 : 110 continue 692 : segdes,descr 693 : irigel(3,nrigel)=descr 694 : * 695 : * creation du segment xmatri 696 : * 697 : nelrig=nelpa2 698 : segini xmatri 699 : irigel(4,nrigel)=xmatri 700 : * 701 : * ce qu'on cree est unilateral 702 : * 703 : irigel(6,nrigel)=1 704 : * 705 : * ce qu'on cree est symetrique 706 : * 707 : irigel(7,nrigel)=0 708 : * 709 : * Nombre d'elements dans la rigidite du type 2 710 : nelri2 = 0 711 : * 712 : DO 120 iel = 1, nelpa2 713 : * 714 : ip1 = lsegm(iel,1) 715 : if (ip1.eq.0) goto 120 716 : C*DBG write (ioimp,*) 'point en relation ligne point ',ip1 717 : lsegm(iel,1) = 0 718 : ip2 = lsegm(iel,2) 719 : jp = lsegm(iel,3) 720 : il1 = lsegm(iel,4) 721 : xn = xsegm(iel,1) 722 : yn = xsegm(iel,2) 723 : zn = xsegm(iel,3) 724 : sn = sqrt(xn*xn+yn*yn+zn*zn) 725 : * 726 : DO 130 jel = iel+1,nelpa2 727 : **pv 728 : goto 130 729 : if (lsegm(jel,1).ne.ip1 .or. lsegm(jel,2).ne.ip2 .or. 730 : & lsegm(jel,3).ne.jp) goto 130 731 : C*DBG write(ioimp,*) 'point en double ' 732 : * verif memes directions 733 : xn1 = xsegm(jel,1) 734 : yn1 = xsegm(jel,2) 735 : zn1 = xsegm(jel,3) 736 : if (xn*xn1+yn*yn1+zn*zn1 .lt. 0.9d0*sn) goto 130 737 : lsegm(jel,1)=0 738 : xn = xn + xn1 739 : yn = yn + yn1 740 : zn = zn + zn1 741 : sn = sqrt(xn*xn+yn*yn+zn*zn) 742 : 130 CONTINUE 743 : if (.not.(sn.lt.1d-10) .and. .not.(sn.gt.1d-10)) 744 : & write(ioimp,*) 'Prob 6.1 - impo32' 745 : if (abs(sn).le.xpetit) write(ioimp,*) 'Prob 6.2 - impo32' 746 : C*DBG write(ioimp,*) ' xn yn zn sn ',xn,yn,zn,sn 747 : xn = xn / sn 748 : yn = yn / sn 749 : zn = zn / sn 750 : if (abs(xn).lt.1d-12) xn = 0.d0 751 : if (abs(yn).lt.1d-12) yn = 0.d0 752 : if (abs(zn).lt.1d-12) zn = 0.d0 753 : * 754 : ipv = (ip1-1)*idimp1 755 : xp1 = xcoor(ipv+1) 756 : yp1 = xcoor(ipv+2) 757 : zp1 = xcoor(ipv+3) 758 : ipv = (ip2-1)*idimp1 759 : xp2 = xcoor(ipv+1) 760 : yp2 = xcoor(ipv+2) 761 : zp2 = xcoor(ipv+3) 762 : ipv = (jp -1)*idimp1 763 : xp = xcoor(ipv+1) 764 : yp = xcoor(ipv+2) 765 : zp = xcoor(ipv+3) 766 : b1 = (xp-xp2)*(xp1-xp2)+(yp-yp2)*(yp1-yp2)+(zp-zp2)*(zp1-zp2) 767 : b2 = (xp-xp1)*(xp2-xp1)+(yp-yp1)*(yp2-yp1)+(zp-zp1)*(zp2-zp1) 768 : b1 = max(b1,0.d0) 769 : b2 = max(b2,0.d0) 770 : bt = b1+b2 771 : if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10)) 772 : & write(ioimp,*) 'Prob 7.1 - impo32' 773 : if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 7.2 - impo32' 774 : b1 = b1 / bt 775 : b2 = b2 / bt 776 : * recalcul de la direction et du jeu en fonction des normales aux sommets 777 : xn = b1 *xms(ip1)+b2 *xms(ip2) 778 : yn = b1 *yms(ip1)+b2 *yms(ip2) 779 : zn = b1 *zms(ip1)+b2 *zms(ip2) 780 : sn = sqrt (xn*xn + yn*yn + zn*zn) 781 : xn=xn/sn 782 : yn=yn/sn 783 : zn=zn/sn 784 : xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn 785 : xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn 786 : xjeu = b1 * xjeu1 + b2 * xjeu2 787 : C*DBG write(ioimp,*) 'type 2 ligne ',b1,b2 788 : * 789 : nelri2 = nelri2 + 1 790 : nelch = nelch + 1 791 : * 792 : * on range dans le meleme 793 : * 794 : num(1,nelri2) = il1 795 : num(2,nelri2) = ip1 796 : num(3,nelri2) = ip2 797 : num(4,nelri2) = jp 798 : * 799 : * on remplit le xmatri 800 : * 801 : * lambda 802 : re( 1,1,nelri2) = 0d0 803 : * ip1 804 : re( 2,1,nelri2) = xn*b1 805 : re( 3,1,nelri2) = yn*b1 806 : re( 4,1,nelri2) = zn*b1 807 : * ip2 808 : re( 5,1,nelri2) = xn*b2 809 : re( 6,1,nelri2) = yn*b2 810 : re( 7,1,nelri2) = zn*b2 811 : * jp 812 : re( 8,1,nelri2) = -xn 813 : re( 9,1,nelri2) = -yn 814 : re(10,1,nelri2) = -zn 815 : * on transpose 816 : do 140 ic=2,10 817 : re(1,ic,nelri2) = re(ic,1,nelri2) 818 : 140 continue 819 : * le reste est nul 820 : * 821 : * remplissage du champoint de depi (jeu) 822 : * 823 : ipt7.num(1,nelch) = il1 824 : vpocha(nelch,1) = xjeu 825 : C*DBG write(ioimp,*) ' jeu ',xjeu,ip1,ip2,jp,il1,il2 826 : * 827 : 120 CONTINUE 828 : * 829 : * Ajustement au plus juste puis desactivation des segments lies 830 : * a la rigidite du type 2 831 : if (nelri2.ne.nelrig) then 832 : nelrig = nelri2 833 : segadj xmatri 834 : nbelem = nelri2 835 : nbnn = 4 836 : segadj,meleme 837 : endif 838 : segdes,xmatri,meleme 839 : * 840 : write (ioimp,*) ' nb relation type 2 ',nelri2 841 : C*DBG if (nelri2.eq.0) irigel(6,nrigel)=0 842 : * 843 : * Relation du type 1 : noeud-noeud 844 : *---------------------------------- 845 : 200 CONTINUE 846 : if (nelpa1.eq.0) goto 300 847 : * 848 : * Il y a au maximum nelpa1 relations de type 1 ! 849 : * 850 : * Ajustement du chpoint de depi (jeu) : mpoval et ipt7 851 : * 852 : n = nelch + nelpa1 853 : segadj,mpoval 854 : nbnn = 1 855 : nbelem = nelch + nelpa1 856 : nbsous = 0 857 : nbref = 0 858 : segadj,ipt7 859 : * 860 : * creation de la raideur de conditions 861 : * 862 : nrigel = nrigel + 1 863 : segadj mrigid 864 : coerig(nrigel) = 1.d0 865 : * 866 : * creation du meleme associe a la relation 867 : * 868 : nbnn = 3 869 : nbelem = nelpa1 870 : nbsous = 0 871 : nbref = 0 872 : segini meleme 873 : itypel = 22 874 : irigel(1,nrigel) = meleme 875 : * 876 : * creation du descriptif commun a toutes les raideurs 877 : * 878 : nligrp=7 879 : nligrd=nligrp 880 : segini descr 881 : lisinc(1)='LX ' 882 : lisdua(1)='FLX ' 883 : noelep(1)=1 884 : noeled(1)=1 885 : do 210 i=2,7,3 886 : lisinc(i )=modepl(1) 887 : lisinc(i+1)=modepl(2) 888 : lisinc(i+2)=modepl(3) 889 : lisdua(i )=moforc(1) 890 : lisdua(i+1)=moforc(2) 891 : lisdua(i+2)=moforc(3) 892 : noelep(i )=(i+5)/3 893 : noelep(i+1)=noelep(i) 894 : noelep(i+2)=noelep(i) 895 : noeled(i )=noelep(i) 896 : noeled(i+1)=noelep(i) 897 : noeled(i+2)=noelep(i) 898 : 210 continue 899 : segdes,descr 900 : irigel(3,nrigel) = descr 901 : * 902 : * creation du segment xmatri 903 : * 904 : nelrig = nelpa1 905 : segini,xmatri 906 : irigel(4,nrigel)=xmatri 907 : * 908 : * ce qu'on cree est unilateral 909 : * 910 : irigel(6,nrigel)=1 911 : * 912 : * ce qu'on cree est symetrique 913 : * 914 : irigel(7,nrigel)=0 915 : * 916 : * Nombre d'elements dans la rigidite du type 1 917 : nelri1 = 0 918 : * 919 : DO 220 iel = 1, nelpa1 920 : * 921 : ip1 = lpoin(iel,1) 922 : if (ip1.eq.0) goto 220 923 : lpoin(iel,1) = 0 924 : jp = lpoin(iel,2) 925 : il1 = lpoin(iel,3) 926 : * 927 : xn = xpoin(iel,1) 928 : yn = xpoin(iel,2) 929 : zn = xpoin(iel,3) 930 : sn = sqrt(xn*xn+yn*yn+zn*zn) 931 : C*DBG write (ioimp,*) 'point en relation point point ', 932 : C*DBG & ip1,jp,il1,xn,yn,zn 933 : DO 230 jel = iel+1, nelpa1 934 : if (lpoin(jel,1).ne.ip1.or.lpoin(jel,2).ne.jp) goto 231 935 : * verif memes directions 936 : xn1 = xpoin(jel,1) 937 : yn1 = xpoin(jel,2) 938 : zn1 = xpoin(jel,3) 939 : **pv if (xn*xn1+yn*yn1+zn*zn1 .lt. 0.9d0*sn) goto 230 940 : C*DBG write (ioimp,*) 'point en double cas 1',ip1,jp 941 : lpoin(jel,1) = 0 942 : xn = xn + xn1 943 : yn = yn + yn1 944 : zn = zn + zn1 945 : sn = sqrt(xn*xn+yn*yn+zn*zn) 946 : 231 continue 947 : **pv 948 : goto 230 949 : if (lpoin(jel,1).ne.jp.or.lpoin(jel,2).ne.ip1) goto 230 950 : * verif memes directions 951 : xn1 = xpoin(jel,1) 952 : yn1 = xpoin(jel,2) 953 : zn1 = xpoin(jel,3) 954 : if (xn*xn1+yn*yn1+zn*zn1 .gt. -0.9d0*sn) goto 230 955 : C*DBG write (ioimp,*) 'point en double cas 2',jp,ip1 956 : lpoin(jel,1) = 0 957 : xn = xn - xn1 958 : yn = yn - yn1 959 : zn = zn - zn1 960 : sn = sqrt(xn*xn+yn*yn+zn*zn) 961 : 230 CONTINUE 962 : C*DBG write (ioimp,*) ' xn yn zn sn ',xn,yn,zn,sn 963 : if (.not.(sn.lt.1d-10).and. .not.(sn.gt.1d-10)) 964 : & write(ioimp,*) 'Prob 8.1 - impo32' 965 : if (abs(sn).le.xpetit) write(ioimp,*) 'Prob 8.2 - impo32' 966 : xn = xn / sn 967 : yn = yn / sn 968 : zn = zn / sn 969 : if (abs(xn).lt.1d-12) xn = 0.d0 970 : if (abs(yn).lt.1d-12) yn = 0.d0 971 : if (abs(zn).lt.1d-12) zn = 0.d0 972 : * 973 : ipv = (ip1-1)*idimp1 974 : xp1 = xcoor(ipv+1) 975 : yp1 = xcoor(ipv+2) 976 : zp1 = xcoor(ipv+3) 977 : ipv = (jp -1)*idimp1 978 : xp = xcoor(ipv+1) 979 : yp = xcoor(ipv+2) 980 : zp = xcoor(ipv+3) 981 : * recalcul de la direction et du jeu en fonction des normales aux sommets 982 : xn = xms(ip1) 983 : yn = yms(ip1) 984 : zn = zms(ip1) 985 : sn = sqrt (xn*xn + yn*yn + zn*zn) 986 : xn=xn/sn 987 : yn=yn/sn 988 : zn=zn/sn 989 : xjeu = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn 990 : * 991 : nelri1 = nelri1 + 1 992 : nelch = nelch + 1 993 : * 994 : * on range dans le meleme 995 : * 996 : num(1,nelri1) = il1 997 : num(2,nelri1) = ip1 998 : num(3,nelri1) = jp 999 : * 1000 : * on remplit le xmatri 1001 : * 1002 : * lambda 1003 : re(1,1,nelri1) = 0d0 1004 : * ip1 1005 : re(2,1,nelri1) = xn 1006 : re(3,1,nelri1) = yn 1007 : re(4,1,nelri1) = zn 1008 : * jp 1009 : re(5,1,nelri1) = -xn 1010 : re(6,1,nelri1) = -yn 1011 : re(7,1,nelri1) = -zn 1012 : * on transpose 1013 : do 240 ic=2,7 1014 : re(1,ic,nelri1) = re(ic,1,nelri1) 1015 : 240 continue 1016 : * le reste est nul 1017 : * 1018 : * remplissage du champoint de depi (jeu) 1019 : ipt7.num(1,nelch) = il1 1020 : vpocha(nelch,1) = xjeu 1021 : C*DBG write(ioimp,*) ' jeu ',xjeu,ip1,jp,il1,il2 1022 : * 1023 : 220 CONTINUE 1024 : * 1025 : * Ajustement au plus juste puis desactivation des segments lies 1026 : * a la rigidite du type 1 1027 : if (nelri1.ne.nelrig) then 1028 : nelrig = nelri1 1029 : segadj xmatri 1030 : nbelem = nelri1 1031 : nbnn = 3 1032 : segadj,meleme 1033 : endif 1034 : segdes,xmatri,meleme 1035 : * 1036 : C*DBG if (nelri1.eq.0) irigel(6,nrigel)=0 1037 : write(ioimp,*) ' nb relation type 1 ',nelri1 1038 : * 1039 : * Fin du traitement de la formulation standard/forte 1040 : * ---------------------------------------------------- 1041 : 300 CONTINUE 1042 : * Destruction des segments locaux devenus inutiles 1043 : segsup,mfopa1,mfopa2,mfopa3 1044 : 1045 : GOTO 1000 1046 : 1047 : *======================================================================= 1048 : *= Formulation "faible" 3D : relation triangle-triangle 1049 : *======================================================================= 1050 : * Element de contact 3D a 7 noeuds (1+6) 1051 : 500 CONTINUE 1052 : * 1053 : * Petit segment de travail 1054 : segini,mfaible 1055 : * 1056 : * Relation du type 0 : triangle-triangle (faible) 1057 : *--------------------------------------- 1058 : * 1059 : * Creation du meleme associe a la relation 1060 : * 1061 : nbnn = 7 1062 : nbelem = incnel 1063 : nbsous = 0 1064 : nbref = 0 1065 : segini,meleme 1066 : itypel = 22 1067 : irigel(1,nrigel) = meleme 1068 : * 1069 : * Creation du descriptif commun a toutes les raideurs 1070 : * 1071 : nligrp = 19 1072 : nligrd = nligrp 1073 : segini,descr 1074 : lisinc(1) = 'LX ' 1075 : lisdua(1) = 'FLX ' 1076 : noelep(1) = 1 1077 : noeled(1) = 1 1078 : do 510 i = 2, nligrp, 3 1079 : lisinc(i ) = modepl(1) 1080 : lisinc(i+1) = modepl(2) 1081 : lisinc(i+2) = modepl(3) 1082 : lisdua(i ) = moforc(1) 1083 : lisdua(i+1) = moforc(2) 1084 : lisdua(i+2) = moforc(3) 1085 : noelep(i ) = (i+5)/3 1086 : noelep(i+1) = noelep(i) 1087 : noelep(i+2) = noelep(i) 1088 : noeled(i ) = noelep(i) 1089 : noeled(i+1) = noelep(i) 1090 : noeled(i+2) = noelep(i) 1091 : 510 continue 1092 : segdes,descr 1093 : irigel(3,nrigel) = descr 1094 : * 1095 : * creation du segment xmatri 1096 : * 1097 : nelrig = incnel 1098 : segini,xmatri 1099 : irigel(4,nrigel) = xmatri 1100 : * 1101 : * ce qu'on cree est unilateral 1102 : * 1103 : irigel(6,nrigel) = 1 1104 : * 1105 : * ce qu'on cree est symetrique 1106 : * 1107 : irigel(7,nrigel) = 0 1108 : * 1109 : * Nombre d'elements crees dans meleme=irigel(1,nrigel), ipt7 et mpoval 1110 : nelri0 = 0 1111 : * 1112 : * Boucle sur les elements du maillage de contact/frottement "faible" 1113 : * 1114 : DO 520 iel = 1, nbel1 1115 : * 1116 : ip1 = ipt1.num(2,iel) 1117 : ip2 = ipt1.num(3,iel) 1118 : ip3 = ipt1.num(4,iel) 1119 : jp1 = ipt1.num(5,iel) 1120 : jp2 = ipt1.num(6,iel) 1121 : jp3 = ipt1.num(7,iel) 1122 : C*DBG write(ioimp,*) iel,ip1,ip2,ip3,jp1,jp2,jp3,ipt1.num(1,iel) 1123 : * 1124 : * Verification (provisoire) qu'il n'y a pas de noeud double 1125 : if (ip1.eq.jp1) goto 520 1126 : if (ip1.eq.jp2) goto 520 1127 : if (ip1.eq.jp3) goto 520 1128 : if (ip2.eq.jp1) goto 520 1129 : if (ip2.eq.jp2) goto 520 1130 : if (ip2.eq.jp3) goto 520 1131 : if (ip3.eq.jp1) goto 520 1132 : if (ip3.eq.jp2) goto 520 1133 : if (ip3.eq.jp3) goto 520 1134 : * 1135 : * Definition du triangle T1(ip1,ip2,ip3) 1136 : * - Coordonnees des noeuds de T1 1137 : ipv = (ip1-1)*idimp1 1138 : xp1 = xcoor(ipv+1) 1139 : yp1 = xcoor(ipv+2) 1140 : zp1 = xcoor(ipv+3) 1141 : ipv = (ip2-1)*idimp1 1142 : xp2 = xcoor(ipv+1) 1143 : yp2 = xcoor(ipv+2) 1144 : zp2 = xcoor(ipv+3) 1145 : ipv = (ip3-1)*idimp1 1146 : xp3 = xcoor(ipv+1) 1147 : yp3 = xcoor(ipv+2) 1148 : zp3 = xcoor(ipv+3) 1149 : * - Barycentre de T1 1150 : xgT1 = (xp1 + xp2 + xp3) * X1s3 1151 : ygT1 = (yp1 + yp2 + yp3) * X1s3 1152 : zgT1 = (zp1 + zp2 + zp3) * X1s3 1153 : * - Normale au triangle T1 1154 : xnT1 = (yp1-yp2)*(zp2-zp3) - (zp1-zp2)*(yp2-yp3) 1155 : ynT1 = (zp1-zp2)*(xp2-xp3) - (xp1-xp2)*(zp2-zp3) 1156 : znT1 = (xp1-xp2)*(yp2-yp3) - (yp1-yp2)*(xp2-xp3) 1157 : dnT1 = SQRT(xnT1*xnT1+ynT1*ynT1+znT1*znT1) 1158 : IF (.not.(dnT1.lt.1d-10) .and. .not.(dnT1.gt.1d-10)) 1159 : & write(ioimp,*) 'FAIBle - 1.1 - impo32' 1160 : IF (abs(dnT1).le.xpetit) write(ioimp,*) 'FAIBle - 1.2 - impo32' 1161 : xnT1 = xnT1 / dnT1 1162 : ynT1 = ynT1 / dnT1 1163 : znT1 = znT1 / dnT1 1164 : if (abs(xnT1).lt.1d-12) xnT1 = 0.d0 1165 : if (abs(ynT1).lt.1d-12) ynT1 = 0.d0 1166 : if (abs(znT1).lt.1d-12) znT1 = 0.d0 1167 : C*DBG write(ioimp,*) 'Normale au plan T1',xnT1,ynT1,znT1,dnT1 1168 : * - Rayon de la sphere "englobante" de l'element (centre=barycentre) 1169 : d1 = (xgT1-xp1)**2 + (ygT1-yp1)**2 + (zgT1-zp1)**2 1170 : d2 = (xgT1-xp2)**2 + (ygT1-yp2)**2 + (zgT1-zp2)**2 1171 : d3 = (xgT1-xp3)**2 + (ygT1-yp3)**2 + (zgT1-zp3)**2 1172 : RayT1 = SQRT(max(d1,d2,d3)) 1173 : * 1174 : * Definition du triangle T2(jp1,jp2,jp3) 1175 : * - Coordonnees des noeuds de T2 1176 : ipv = (jp1-1)*idimp1 1177 : xq1 = xcoor(ipv+1) 1178 : yq1 = xcoor(ipv+2) 1179 : zq1 = xcoor(ipv+3) 1180 : ipv = (jp2-1)*idimp1 1181 : xq2 = xcoor(ipv+1) 1182 : yq2 = xcoor(ipv+2) 1183 : zq2 = xcoor(ipv+3) 1184 : ipv = (jp3-1)*idimp1 1185 : xq3 = xcoor(ipv+1) 1186 : yq3 = xcoor(ipv+2) 1187 : zq3 = xcoor(ipv+3) 1188 : * - Barycentre de T2 1189 : xgT2 = (xq1 + xq2 + xq3) * X1s3 1190 : ygT2 = (yq1 + yq2 + yq3) * X1s3 1191 : zgT2 = (zq1 + zq2 + zq3) * X1s3 1192 : * - Normale au triangle T2 1193 : xnT2 = (yq1-yq2)*(zq2-zq3) - (zq1-zq2)*(yq2-yq3) 1194 : ynT2 = (zq1-zq2)*(xq2-xq3) - (xq1-xq2)*(zq2-zq3) 1195 : znT2 = (xq1-xq2)*(yq2-yq3) - (yq1-yq2)*(xq2-xq3) 1196 : dnT2 = SQRT(xnT2*xnT2+ynT2*ynT2+znT2*znT2) 1197 : IF (.not.(dnT2.lt.1d-10) .and. .not.(dnT2.gt.1d-10)) 1198 : & write(ioimp,*) 'FAIBle - 2.1 - impo32' 1199 : IF (abs(dnT2).le.xpetit) write(ioimp,*) 'FAIBle - 2.2 - impo32' 1200 : xnT2 = xnT2 / dnT2 1201 : ynT2 = ynT2 / dnT2 1202 : znT2 = znT2 / dnT2 1203 : if (abs(xnT2).lt.1d-12) xnT2 = 0.d0 1204 : if (abs(ynT2).lt.1d-12) ynT2 = 0.d0 1205 : if (abs(znT2).lt.1d-12) znT2 = 0.d0 1206 : C*DBG write(ioimp,*) 'Normale au plan T2',xnT2,ynT2,znT2,dnT2 1207 : * - Rayon de la sphere "englobante" de l'element (centre=barycentre) 1208 : d1 = (xgT2-xq1)**2 + (ygT2-yq1)**2 + (zgT2-zq1)**2 1209 : d2 = (xgT2-xq2)**2 + (ygT2-yq2)**2 + (zgT2-zq2)**2 1210 : d3 = (xgT2-xq3)**2 + (ygT2-yq3)**2 + (zgT2-zq3)**2 1211 : RayT2 = SQRT(max(d1,d2,d3)) 1212 : 1213 : * Orientation respective correcte des 2 normales 1214 : * Les triangles T1 et T2 etant decrits dans le meme sens ("prisme"), 1215 : * le produit scalaire de leur normale doit etre positif ou nul ! 1216 : scal = xnT1 * xnT2 + ynT1 * ynT2 + znT1 * znT2 1217 : C*DBG write(ioimp,*) 'Prod.scal des normales',scal 1218 : if (scal.lt.0.) goto 520 1219 : * 1220 : * Distance entre les centres de gravite 1221 : * Les spheres "englobantes" doivent s'intersecter pour contact potentiel 1222 : dist = SQRT((xgT2-xgT1)**2 + (ygT2-ygT1)**2 + (zgT2-zgT1)**2) 1223 : IF (dist.GT.RayT1+RayT2) goto 520 1224 : * 1225 : * Critere d'acceptation de l'element de contact 1226 : * 1227 : * Definition du plan de contact (surface intermediaire) 1228 : * - Normale (unitaire) du plan de contact 1229 : xnC = xnT1 + xnT2 1230 : ynC = ynT1 + ynT2 1231 : znC = znT1 + znT2 1232 : dnC = SQRT(xnC*xnC + ynC*ynC + znC*znC) 1233 : IF (.not.(dnC.lt.1d-10) .and. .not.(dnC.GT.1D-10)) 1234 : & write(ioimp,*) 'FAIBle - 3.1 - impo32' 1235 : IF (abs(dnC).le.xpetit) write(ioimp,*) 'FAIBle - 3.2 - impo32' 1236 : xnC = xnC / dnC 1237 : ynC = ynC / dnC 1238 : znC = znC / dnC 1239 : if (abs(xnC).lt.1d-12) xnC = 0.d0 1240 : if (abs(ynC).lt.1d-12) ynC = 0.d0 1241 : if (abs(znC).lt.1d-12) znC = 0.d0 1242 : C*DBG write(ioimp,*) 'Normale au plan Contact',xnC,ynC,znC,dnC 1243 : 1244 : * - Point definissant le plan de contact 1245 : * = milieu des barycentres des triangles T1 et T2 1246 : xgC = (xgT1 + xgT2) * 0.5 1247 : ygC = (ygT1 + ygT2) * 0.5 1248 : zgC = (zgT1 + zgT2) * 0.5 1249 : * - "Distance" de reference au plan de contact 1250 : dgC = xgC*xnC + ygC*ynC + zgC*znC 1251 : * 1252 : * Distance "signee" des triangles T1 et T2 au plan de contact 1253 : dp1C = (xp1*xnC + yp1*ynC + zp1*znC) - dgC 1254 : dp2C = (xp2*xnC + yp2*ynC + zp2*znC) - dgC 1255 : dp3C = (xp3*xnC + yp3*ynC + zp3*znC) - dgC 1256 : dq1C = (xq1*xnC + yq1*ynC + zq1*znC) - dgC 1257 : dq2C = (xq2*xnC + yq2*ynC + zq2*znC) - dgC 1258 : dq3C = (xq3*xnC + yq3*ynC + zq3*znC) - dgC 1259 : * Projection des triangles T1 et T2 sur le plan de contact 1260 : * T1C = triangle T1 projete sur plan de Contact 1261 : cPT1C(1,1) = xp1 - dp1C * xnC 1262 : cPT1C(1,2) = yp1 - dp1C * ynC 1263 : cPT1C(1,3) = zp1 - dp1C * znC 1264 : cPT1C(2,1) = xp2 - dp2C * xnC 1265 : cPT1C(2,2) = yp2 - dp2C * ynC 1266 : cPT1C(2,3) = zp2 - dp2C * znC 1267 : cPT1C(3,1) = xp3 - dp3C * xnC 1268 : cPT1C(3,2) = yp3 - dp3C * ynC 1269 : cPT1C(3,3) = zp3 - dp3C * znC 1270 : * T2C = triangle T2 projete sur plan de Contact 1271 : cPT2C(1,1) = xq1 - dq1C * xnC 1272 : cPT2C(1,2) = yq1 - dq1C * ynC 1273 : cPT2C(1,3) = zq1 - dq1C * znC 1274 : cPT2C(2,1) = xq2 - dq2C * xnC 1275 : cPT2C(2,2) = yq2 - dq2C * ynC 1276 : cPT2C(2,3) = zq2 - dq2C * znC 1277 : cPT2C(3,1) = xq3 - dq3C * xnC 1278 : cPT2C(3,2) = yq3 - dq3C * ynC 1279 : cPT2C(3,3) = zq3 - dq3C * znC 1280 : * 1281 : * On determine quelle est la composante maximale de la normale pour 1282 : * projeter les triangles T1C et T2C sur un plan parallele aux axes 1283 : * ("le plus orthogonal" a cette normale), ce qui maximise leur aire. 1284 : r_x = abs(xnC) 1285 : r_y = abs(ynC) 1286 : r_z = abs(znC) 1287 : if (r_x.gt.r_y) then 1288 : inC = 1 1289 : if (r_x.lt.r_z) inC = 3 1290 : else 1291 : inC = 2 1292 : if (r_y.lt.r_z) inC = 3 1293 : endif 1294 : * Projection des triangles T1C et T2C sur ce plan // aux axes (A) 1295 : * T1A = triangle T1C projete sur ce plan (A) 1296 : * T2A = triangle T2C projete sur ce plan (A) 1297 : * On prend soin d'orienter les triangles T1A et T2A dans le sens 1298 : * direct (en tenant compte du signe de la composante normale !) 1299 : GOTO (531,532,533), inC 1299 : call erreur(5) 1300 : return 1301 : * Normale NC selon "x" : A =(y,z) 1302 : 531 if (xnC.lt.0.) inC = -inC 1303 : inX = 2 1304 : inY = 3 1305 : goto 534 1306 : * Normale NC selon "y" : A =(x,z) 1307 : 532 if (ynC.lt.0.) inC = -inC 1308 : inX = 1 1309 : inY = 3 1310 : goto 534 1311 : * Normale NC selon "z" : A =(x,y) 1312 : 533 if (znC.lt.0.) inC = -inC 1313 : inX = 1 1314 : inY = 2 1315 : goto 534 1316 : 534 continue 1317 : cPT1A(1,1) = cPT1C(1,inX) 1318 : cPT1A(1,2) = cPT1C(1,inY) 1319 : cPT2A(1,1) = cPT2C(1,inX) 1320 : cPT2A(1,2) = cPT2C(1,inY) 1321 : if (inC.gt.0) then 1322 : cPT1A(2,1) = cPT1C(2,inX) 1323 : cPT1A(2,2) = cPT1C(2,inY) 1324 : cPT2A(2,1) = cPT2C(2,inX) 1325 : cPT2A(2,2) = cPT2C(2,inY) 1326 : cPT1A(3,1) = cPT1C(3,inX) 1327 : cPT1A(3,2) = cPT1C(3,inY) 1328 : cPT2A(3,1) = cPT2C(3,inX) 1329 : cPT2A(3,2) = cPT2C(3,inY) 1330 : else 1331 : cPT1A(2,1) = cPT1C(3,inX) 1332 : cPT1A(2,2) = cPT1C(3,inY) 1333 : cPT2A(2,1) = cPT2C(3,inX) 1334 : cPT2A(2,2) = cPT2C(3,inY) 1335 : cPT1A(3,1) = cPT1C(2,inX) 1336 : cPT1A(3,2) = cPT1C(2,inY) 1337 : cPT2A(3,1) = cPT2C(2,inX) 1338 : cPT2A(3,2) = cPT2C(2,inY) 1339 : endif 1340 : cPT1A(4,1) = cPT1A(1,1) 1341 : cPT1A(4,2) = cPT1A(1,2) 1342 : * Calcul des cotes des triangles T1A et T2A (12,23,31) 1343 : * Utile pour connaitre ensuite la normale aux cotes des triangles 1344 : vT1A(1,1) = cPT1A(2,1) - cPT1A(1,1) 1345 : vT1A(1,2) = cPT1A(2,2) - cPT1A(1,2) 1346 : vT1A(2,1) = cPT1A(3,1) - cPT1A(2,1) 1347 : vT1A(2,2) = cPT1A(3,2) - cPT1A(2,2) 1348 : vT1A(3,1) = cPT1A(1,1) - cPT1A(3,1) 1349 : vT1A(3,2) = cPT1A(1,2) - cPT1A(3,2) 1350 : C*nu vT2A(1,1) = cPT2A(2,1) - cPT2A(1,1) 1351 : vT2A(1,2) = cPT2A(2,2) - cPT2A(1,2) 1352 : C*nu vT2A(2,1) = cPT2A(3,1) - cPT2A(2,1) 1353 : vT2A(2,2) = cPT2A(3,2) - cPT2A(2,2) 1354 : C*nu vT2A(3,1) = cPT2A(1,1) - cPT2A(3,1) 1355 : vT2A(3,2) = cPT2A(1,2) - cPT2A(3,2) 1356 : * Calcul de la surface de chacun des triangles T1A et T2A 1357 : * (en fait, on calcule le double, mais par la suite on ne considere que 1358 : * des rapports de surfaces de triangle) 1359 : SuT1A = (cPT1A(2,1)+cPT1A(1,1)) * vT1A(1,2) 1360 : & + (cPT1A(3,1)+cPT1A(2,1)) * vT1A(2,2) 1361 : & + (cPT1A(1,1)+cPT1A(3,1)) * vT1A(3,2) 1362 : SuT2A = (cPT2A(2,1)+cPT2A(1,1)) * vT2A(1,2) 1363 : & + (cPT2A(3,1)+cPT2A(2,1)) * vT2A(2,2) 1364 : & + (cPT2A(1,1)+cPT2A(3,1)) * vT2A(3,2) 1365 : C*DBG write(ioimp,*) 'Surfaces T1A - T2A :', 0.5*SuT1A,0.5*SuT2A 1366 : C*DBG if (SuT1A.gt.100.*SuT2A .or. SuT2A.gt.100.*SuT1A) 1367 : C*DBG & write(ioimp,*) 'Rapport des surfaces tres important !' 1368 : 1369 : * On initialise l'intersection des triangles avec T2A 1370 : nPIn0 = 3 1371 : do 540 i = 1, nPIn0 1372 : cPIn0(i,1) = cPT2A(i,1) 1373 : cPIn0(i,2) = cPT2A(i,2) 1374 : 540 continue 1375 : * Determination de l'intersection des triangles T1A et T2A 1376 : * en regardant progressivement l'intersection avec les cotes de T1A 1377 : * Note : On utilise le fait que les polygones traites sont convexes ! 1378 : DO 550 i = 1, 3 1379 : vN1x = -vT1A(i,2) 1380 : vN1y = vT1A(i,1) 1381 : scal = vN1x * cPT1A(i,1) + vN1y * cPT1A(i,2) 1382 : * 1383 : ipos = 0 1384 : ineg = 0 1385 : jpos = 0 1386 : jneg = 0 1387 : do 551 j = 1, nPIn0 1388 : test(j) = (vN1x * cPIn0(j,1) + vN1y * cPIn0(j,2)) - scal 1389 : if (test(j).gt.1.d-10) then 1390 : * if (test(j).gt.0.) then 1391 : ipos = ipos + 1 1392 : if (jneg.ne.0 .and. jpos.eq.0) jpos = j 1393 : else if (test(j).lt.-1.d-10) then 1394 : * else if (test(j).lt.0.) then 1395 : ineg = ineg + 1 1396 : if (jneg.eq.0) jneg = j 1397 : else 1398 : test(j) = 0. 1399 : endif 1400 : 551 continue 1401 : * 1402 : * 1) Cas ou ipos = 0 : 1403 : * 1.1) il n'y pas d'intersection (ineg=nPIn0) 1404 : * 1.2) l'intersection se limite a 1 point (ineg=nPIn0-1) 1405 : * 1.3) l'intersection a 1 segment sur une frontiere (ineg=nPIn0-2) 1406 : * La surface correspondante est nulle et il n'y a donc pas de contact ! 1407 : if (ipos.eq.0) then 1408 : C*DBG write(ioimp,*) 'Intersection a ',nPIn0-ineg,' points' 1409 : goto 520 1410 : endif 1411 : * 2) Cas ou ipos > 0 : Il y a une intersection a calculer 1412 : * 2.1) ineg = 0 : Tous les points sont "du bon cote" et l'intersection 1413 : * correspond au polygone de depart ! 1414 : if (ineg.eq.0) then 1415 : C*DBG write(ioimp,*) 'Intersection du bon cote - rien a faire' 1416 : goto 550 1417 : endif 1418 : * 2.2) ineg > 0 : Il y a deux intersections a determiner car il y a 1419 : * au moins un point du "mauvais cote" 1420 : jpos = max(1,jpos) 1421 : jpr = jpos - 1 1422 : if (jpos.eq.1) jpr = nPIn0 1423 : nPIn = 1 1424 : if (test(jpr).lt.0.) then 1425 : r_z = test(jpos) / (test(jpos) - test(jpr)) 1426 : cPIn(nPIn,1) = cPIn0(jpos,1) 1427 : & + r_z*(cPIn0(jpr,1) - cPIn0(jpos,1)) 1428 : cPIn(nPIn,2) = cPIn0(jpos,2) 1429 : & + r_z*(cPIn0(jpr,2) - cPIn0(jpos,2)) 1430 : else 1431 : cPIn(nPIn,1) = cPIn0(jpr,1) 1432 : cPIn(nPIn,2) = cPIn0(jpr,2) 1433 : endif 1434 : do 552 j = 1, ipos 1435 : nPIn = nPIn + 1 1436 : cPIn(nPIn,1) = cPIn0(jpos,1) 1437 : cPIn(nPIn,2) = cPIn0(jpos,2) 1438 : jpr = jpos 1439 : jpos = jpos + 1 1440 : if (jpos.gt.nPIn0) jpos = 1 1441 : 552 continue 1442 : nPIn = nPIn + 1 1443 : if (test(jpos).lt.0.) then 1444 : r_z = test(jpr) / (test(jpr) - test(jpos)) 1445 : cPIn(nPIn,1) = cPIn0(jpr,1) 1446 : & + r_z*(cPIn0(jpos,1) - cPIn0(jpr,1)) 1447 : cPIn(nPIn,2) = cPIn0(jpr,2) 1448 : & + r_z*(cPIn0(jpos,2) - cPIn0(jpr,2)) 1449 : else 1450 : cPIn(nPIn,1) = cPIn0(jpos,1) 1451 : cPIn(nPIn,2) = cPIn0(jpos,2) 1452 : endif 1453 : * Mise a jour cPIn0 pour la suite du traitement des intersections... 1454 : nPIn0 = nPIn 1455 : do 554 j = 1, nPIn0 1456 : cPIn0(j,1) = cPIn(j,1) 1457 : cPIn0(j,2) = cPIn(j,2) 1458 : 554 continue 1459 : C*DBG-F write(ioimp,*) 'Intersection ',i,' a ',nPIn0,' points' 1460 : C*DBG-F do j = 1, nPIn0 1461 : C*DBG-F write(ioimp,*) ' ',j,nbpts0+6+j,cPIn0(j,1),cPIn0(j,2) 1462 : C*DBG-F enddo 1463 : * 1464 : 550 CONTINUE 1465 : * 1466 : * Calcul de la surface de l'intersection et de son centre de gravite 1467 : r_z = cPIn0(nPIn0,1)*cPIn0(1,2) - cPIn0(1,1)*cPIn0(nPIn0,2) 1468 : SuPIn = r_z 1469 : xGIn = (cPIn0(nPIn0,1)+cPIn0(1,1)) * r_z 1470 : yGIn = (cPIn0(nPIn0,2)+cPIn0(1,2)) * r_z 1471 : do 560 i = 1, nPIn0-1 1472 : j = i + 1 1473 : r_z = cPIn0(i,1)*cPIn0(j,2) - cPIn0(j,1)*cPIn0(i,2) 1474 : SuPIn = SuPIn + r_z 1475 : xGIn = xGIn + (cPIn0(i,1)+cPIn0(j,1)) * r_z 1476 : yGIn = yGIn + (cPIn0(i,2)+cPIn0(j,2)) * r_z 1477 : 560 continue 1478 : if (SuPIn .lt. 1.E-4*max(SuT1A,SuT2A) ) then 1479 : C*DBG write(ioimp,*) 'Intersection a une surface negligeable !' 1480 : goto 520 1481 : endif 1482 : r_z = X1s3 / SupIn 1483 : SuPIn = 0.5 * SupIn 1484 : xGIn = xGIn * r_z 1485 : yGIn = yGIn * r_z 1486 : C*DBG write(ioimp,*) 'Surface Intersection :',SuPIn 1487 : * 1488 : * Calcul des coordonnees barycentriques du centre de gravite 1489 : * pour le triangle T1A 1490 : xpip1 = xGIn - cPT1A(1,1) 1491 : ypip1 = yGIn - cPT1A(1,2) 1492 : xpip2 = xGIn - cPT1A(2,1) 1493 : ypip2 = yGIn - cPT1A(2,2) 1494 : xpip3 = xGIn - cPT1A(3,1) 1495 : ypip3 = yGIn - cPT1A(3,2) 1496 : b1T1 = xpip2 * ypip3 - ypip2 * xpip3 1497 : b2T1 = xpip3 * ypip1 - ypip3 * xpip1 1498 : b3T1 = xpip1 * ypip2 - ypip1 * xpip2 1499 : bt = b1T1 + b2T1 + b3T1 1500 : if (abs(bt-SuT1A) .gt. 1.E-3) 1501 : & write(ioimp,*) 'Prob bt-SuT1A',bt,SuT1A 1502 : b1T1 = b1T1 / SuT1A 1503 : b2T1 = b2T1 / SuT1A 1504 : b3T1 = b3T1 / SuT1A 1505 : * 1506 : * pour le triangle T2A 1507 : xpip1 = xGIn - cPT2A(1,1) 1508 : ypip1 = yGIn - cPT2A(1,2) 1509 : xpip2 = xGIn - cPT2A(2,1) 1510 : ypip2 = yGIn - cPT2A(2,2) 1511 : xpip3 = xGIn - cPT2A(3,1) 1512 : ypip3 = yGIn - cPT2A(3,2) 1513 : b1T2 = xpip2 * ypip3 - ypip2 * xpip3 1514 : b2T2 = xpip3 * ypip1 - ypip3 * xpip1 1515 : b3T2 = xpip1 * ypip2 - ypip1 * xpip2 1516 : bt = b1T2 + b2T2 + b3T2 1517 : if (abs(bt-SuT2A) .gt. 1.E-3) 1518 : & write(ioimp,*) 'Prob bt-SuT2A',bt,SuT2A 1519 : b1T2 = b1T2 / SuT2A 1520 : b2T2 = b2T2 / SuT2A 1521 : b3T2 = b3T2 / SuT2A 1522 : * 1523 : * Calcul du jeu 1524 : xjeu = (b1T2*dq1C + b2T2*dq2C + b3T2*dq3C) 1525 : & - (b1T1*dp1C + b2T1*dp2C + b3T1*dp3C) 1526 : C*DBG write(ioimp,*) 'Jeu =',xjeu 1527 : * 1528 : * Ajout d'une relation triangle-triangle 1529 : * 1530 : nelri0 = nelri0 + 1 1531 : nelch = nelch + 1 1532 : * 1533 : * on ajuste les differents segments si necesaire 1534 : * 1535 : if (nelri0.gt.nelrig) then 1536 : nelrig = nelrig + incnel 1537 : segadj,xmatri 1538 : nbelem = nbelem + incnel 1539 : nbnn = 7 1540 : segadj,meleme 1541 : nbnn = 1 1542 : segadj,ipt7 1543 : n = n + incnel 1544 : segadj,mpoval 1545 : endif 1546 : * 1547 : * Mise a jour du meleme 1548 : * 1549 : if (ipt1.num(1,iel).eq.ILGNI) then 1550 : segact ipt1*mod 1551 : nbpts=xcoor(/1)/(idim+1)+1 1552 : segadj mcoord 1553 : ipt1.num(1,iel)=nbpts 1554 : endif 1555 : num(1,nelri0) = ipt1.num(1,iel) 1556 : num(2,nelri0) = ip1 1557 : num(3,nelri0) = ip2 1558 : num(4,nelri0) = ip3 1559 : num(5,nelri0) = jp1 1560 : num(6,nelri0) = jp2 1561 : num(7,nelri0) = jp3 1562 : *$ 1563 : * Mise a jour de xmatri 1564 : * 1565 : * lambda 1566 : re( 1,1,nelri0) = 0.d0 1567 : * ip1 1568 : re( 2,1,nelri0) = xnC * b1T1 1569 : re( 3,1,nelri0) = ynC * b1T1 1570 : re( 4,1,nelri0) = znC * b1T1 1571 : * ip2 1572 : re( 5,1,nelri0) = xnC * b2T1 1573 : re( 6,1,nelri0) = ynC * b2T1 1574 : re( 7,1,nelri0) = znC * b2T1 1575 : * ip3 1576 : re( 8,1,nelri0) = xnC * b3T1 1577 : re( 9,1,nelri0) = ynC * b3T1 1578 : re(10,1,nelri0) = znC * b3T1 1579 : * jp1 1580 : re(11,1,nelri0) = -xnC * b1T2 1581 : re(12,1,nelri0) = -ynC * b1T2 1582 : re(13,1,nelri0) = -znC * b1T2 1583 : * jp2 1584 : re(14,1,nelri0) = -xnC * b2T2 1585 : re(15,1,nelri0) = -ynC * b2T2 1586 : re(16,1,nelri0) = -znC * b2T2 1587 : * jp3 1588 : re(17,1,nelri0) = -xnC * b3T2 1589 : re(18,1,nelri0) = -ynC * b3T2 1590 : re(19,1,nelri0) = -znC * b3T2 1591 : * on transpose 1592 : do 580 ic = 2, nligrp 1593 : re(1,ic,nelri0) = re(ic,1,nelri0) 1594 : 580 continue 1595 : * le reste est nul 1596 : * 1597 : * remplissage du champoint de depi (jeu) 1598 : * 1599 : ipt7.num(1,nelch) = ipt1.num(1,iel) 1600 : vpocha(nelch,1) = xjeu 1600 : C*DBG-F call prfaible(ipt1,iel,nPIn0,mfaible) 1601 : * 1602 : 520 CONTINUE 1603 : * 1604 : * Ajustement au plus juste puis desactivation des segments lies 1605 : * a la rigidite du type 0 1606 : if (nelri0.ne.nelrig) then 1607 : nelrig = nelri0 1608 : segadj,xmatri 1609 : nbelem = nelri0 1610 : nbnn = 7 1611 : segadj,meleme 1612 : endif 1613 : segdes,meleme,xmatri 1614 : * 1615 : write(ioimp,*) ' nb relation type 0 ',nelri0 1616 : * if (nelri0.eq.0) irigel(6,nrigel)=0 1617 : * 1618 : * Destruction des segments locaux devenus inutiles 1619 : segsup,mfaible 1620 : * 1621 : C* GOTO 1000 1622 : 1623 : *======================================================================= 1624 : *= Fin du sous-programme : 1625 : *======================================================================= 1626 : 1000 CONTINUE 1627 : * 1628 : * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7 1629 : * puis desactivation du chpoint 1630 : if (n.ne.nelch) then 1631 : n = nelch 1632 : segadj,mpoval 1633 : nbnn = 1 1634 : nbelem = nelch 1635 : nbsous = 0 1636 : nbref = 0 1637 : segadj,ipt7 1638 : endif 1639 : segdes,mpoval,ipt7 1640 : segdes,msoupo,mchpoi 1641 : * Desctivation de la matrice de raideur de contact 1642 : segdes,mrigid 1643 : * 1644 : * Ecriture des resultats de l'operateur : 1644 : call ecrobj('RIGIDITE',mrigid) 1644 : call ecrobj('CHPOINT ',mchpoi) 1645 : * 1646 : 900 CONTINUE 1647 : segdes,ipt1 1648 : 1649 : return 1650 : end 1651 : 1652 : 1653 : 1654 :
© Cast3M 2003 - All rights reserved.
Disclaimer