1 : C IMPOS7 SOURCE CHAT 09/10/09 21:19:02 6519 2 : subroutine impos7 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : -INC CCREEL 6 : -INC CCOPTIO 7 : -INC SMCOORD 8 : -INC SMELEME 9 : -INC SMRIGID 10 : -INC SMCHPOI 11 : dimension vec(3) 12 : character*4 modepl(4),moforc(4),mnez(3),mcle(3) 13 : data modepl /'UX ','UY ','UR ','UZ '/ 14 : data moforc /'FX ','FY ','FR ','FZ '/ 15 : data mnez /'PLAT','HEMI','CONE'/ 16 : data mcle /'MAIT','VECT','ANGL'/ 17 : * 18 : * 19 : * 20 : * lecture des options 21 : * 21 : call lirmot(mcle,3,ire,0) 22 : * 23 : * creation de l'objet MAILLAGE contenant tous les contacts possibles 24 : * 25 : if (ire.eq.1) then 25 : call impos8 26 : return 27 : endif 28 : * 29 : * lecture du maillage support des conditions 30 : * 30 : call lirobj('MAILLAGE',ipt1,1,iretou) 31 : if (ierr.ne.0) return 32 : segact ipt1 32 : if (ipt1.itypel.ne.22) call erreur(851) 33 : * 34 : * cas du contact ponctuel 35 : * 36 : if(ipt1.num(/1).eq.3) then 37 : * 38 : * on test que le dernier noeud est bien identique 39 : * 40 : jp=ipt1.num(3,1) 41 : do 80 iel=2,ipt1.num(/2) 41 : if(ipt1.num(3,iel).ne.jp) call erreur(851) 42 : 80 continue 42 : call lirmot(mnez,3,iret,0) 43 : * 44 : * projectile sans nez 45 : * 46 : if (iret.eq.0) then 46 : call lirobj('POINT',noeud,1,iretou) 47 : do 1 ia=1,idim 48 : vec(ia)=xcoor((noeud-1)*(idim+1)+ia) 49 : 1 CONTINUE 50 : endif 51 : * 52 : * projectile a nez plat 53 : * 54 : if (iret.eq.1) then 54 : call lirree(vlarg,0,iretou) 54 : call lirmot(mcle,3,ire,0) 55 : if(ire.eq.2) then 55 : call lirobj('POINT',noeud,0,iretou) 56 : do 2 ia=1,idim 57 : vec(ia)=xcoor((noeud-1)*(idim+1)+ia) 58 : 2 CONTINUE 59 : else 59 : call erreur(20) 60 : endif 61 : endif 62 : * 63 : * projectile a nez hemispherique 64 : * 65 : if (iret.eq.2) then 65 : call lirree(vrayo,0,iretou) 65 : call lirmot(mcle,3,ire,0) 66 : if(ire.eq.2) then 66 : call lirobj('POINT',noeud,0,iretou) 67 : do 3 ia=1,idim 68 : vec(ia)=xcoor((noeud-1)*(idim+1)+ia) 69 : 3 CONTINUE 70 : else 70 : call erreur(20) 71 : endif 72 : endif 73 : * 74 : * projectile a nez conique 75 : * 76 : if (iret.eq.3) then 76 : call lirree(vlarg,0,iretou) 76 : call lirmot(mcle,3,ire,0) 77 : if(ire.eq.3) then 77 : call lirree(vangl,0,iretou) 78 : vangl=(XPI/2.D0)-(XPI*vangl/180) 79 : else 79 : call erreur(15) 80 : endif 80 : call lirmot(mcle,3,ire,0) 81 : if(ire.eq.2) then 81 : call lirobj('POINT',noeud,0,iretou) 82 : do 4 ia=1,idim 83 : vec(ia)=xcoor((noeud-1)*(idim+1)+ia) 84 : 4 CONTINUE 85 : else 85 : call erreur(20) 86 : endif 87 : endif 88 : endif 89 : * 90 : * creation de la raideur de conditions 91 : * 92 : nrige=7 93 : nrigel=1 94 : segini mrigid 95 : ichole=0 96 : imgeo1=0 97 : imgeo2=0 98 : isupeq=0 99 : iforig=ifomod 100 : coerig(1)=1.d0 101 : mtymat='RIGIDITE' 102 : * 103 : * creation du meleme associe a la relation 104 : * 105 : nbsous=0 106 : nbref=0 107 : nbnn=4 108 : nbelem=0 109 : segini meleme 110 : itypel=22 111 : irigel(1,1)=meleme 112 : * 113 : * impact ligne-ligne (vrai) ou ligne-noeud (sinon) ? 114 : * 115 : if(ipt1.num(/1).eq.4) then 116 : * 117 : * creation des descriptifs communs a toutes les raideurs 118 : * 119 : imo=1 120 : if(ifour.eq.0) imo=3 121 : nligrp=7 122 : nligrd=nligrp 123 : segini descr 124 : irigel(3,1)=descr 125 : lisinc(1)='LX' 126 : lisdua(1)='FLX' 127 : noelep(1)=1 128 : noeled(1)=1 129 : do 5 i=2,7,2 130 : lisinc(i)=modepl(imo) 131 : lisinc(i+1)=modepl(imo+1) 132 : lisdua(i)=moforc(imo) 133 : lisdua(i+1)=moforc(imo+1) 134 : noelep(i)=(i+2)/2 135 : noelep(i+1)=noelep(i) 136 : noeled(i)=noelep(i) 137 : noeled(i+1)=noelep(i) 138 : 5 continue 139 : * 140 : * creation du segment imatri 141 : * 142 : nelrig=50 143 : nelri=0 144 : segini xmatri 145 : irigel(4,1)=xmatri 146 : * 147 : * ce qu'on cree est une égalité 148 : * 149 : irigel(6,1)=0 150 : * 151 : * ce qu'on cree est symetrique 152 : * 153 : irigel(7,1)=0 154 : * 155 : * boucle sur le maillage support 156 : do 10 iel=1,ipt1.num(/2) 157 : * 158 : * noeud maitre 159 : * 160 : ip1=ipt1.num(2,iel) 161 : xp1=xcoor((ip1-1)*(idim+1)+1) 162 : yp1=xcoor((ip1-1)*(idim+1)+2) 163 : * 164 : * deuxieme noeud maitre 165 : * 166 : ip2=ipt1.num(3,iel) 167 : xp2=xcoor((ip2-1)*(idim+1)+1) 168 : yp2=xcoor((ip2-1)*(idim+1)+2) 169 : xr=xp2-xp1 170 : yr=yp2-yp1 171 : * 172 : * noeud esclave 173 : * 174 : jp=ipt1.num(4,iel) 175 : xp=xcoor((jp-1)*(idim+1)+1) 176 : yp=xcoor((jp-1)*(idim+1)+2) 177 : * 178 : * test de recherche du contact 179 : * 180 : if(xp2.gt.xp1) then 181 : if(xp.lt.xp1 .or. xp.gt.xp2) goto 10 182 : else 183 : if(xp.lt.xp2 .or. xp.gt.xp1) goto 10 184 : endif 185 : if((xp-xp1)*yr-(yp-yp1)*xr.gt.0.d0) goto 10 186 : * 187 : * determination du noeud defenseur fictif 188 : * 189 : sr2=xr**2+yr**2 190 : xa=xp2-xp 191 : ya=yp2-yp 192 : alpha=(xa*xr+ya*yr)/sr2 193 : if(abs(alpha).lt.1D-3) alpha=0.D0 194 : if(alpha.gt.0.999D0) alpha=1.D0 195 : if(alpha.lt.-0.999D0) alpha=-1.D0 196 : sr=sqrt(sr2) 197 : xr=xr/sr 198 : yr=yr/sr 199 : * 200 : * creation du xmatri 201 : * 202 : * segini xmatri 203 : * imatri=irigel(4,1) 204 : nelri=nelri+1 205 : if(nelri.gt.nelrig) then 206 : nelrig=nelrig+50 207 : segadj xmatri 208 : endif 209 : * imattt(nelrig)=xmatri 210 : * 211 : * multiplicateur 212 : * 213 : re(1,1,nelri)= 0.d0 214 : * 215 : * 1er noeud maitre ip1 216 : * 217 : re(2,1,nelri)= yr * alpha 218 : re(3,1,nelri)= -xr * alpha 219 : * 220 : * 2eme noeud maitre ip2 221 : * 222 : re(4,1,nelri)= yr * (1-alpha) 223 : re(5,1,nelri)= -xr * (1-alpha) 224 : * 225 : * noeud esclave jp 226 : re(6,1,nelri)= -yr 227 : re(7,1,nelri)= xr 228 : * 229 : * test de conditionnement 230 : do 90 il=1,7 231 : if(abs(re(il,1,nelri)).lt.1D-6) re(il,1,nelri)=0.D0 232 : if(re(il,1,nelri).gt.(1.D0-1D-6)) re(il,1,nelri)=1.D0 233 : if(re(il,1,nelri).lt.(-1.D0+1D-6)) re(il,1,nelri)=-1.D0 234 : 90 continue 235 : * 236 : * transposition du xmatri C 237 : do 30 il=1,1 238 : do 301 ic=1,7 239 : re(il,ic,nelri)=re(ic,il,nelri) 240 : 301 continue 241 : 30 continue 242 : * 243 : * modification du meleme 244 : * 245 : nbelem=num(/2)+1 246 : nbnn=4 247 : nbsous=0 248 : nbref=0 249 : segadj meleme 250 : num(1,nbelem)=ipt1.num(1,iel) 251 : num(2,nbelem)=ip1 252 : num(3,nbelem)=ip2 253 : num(4,nbelem)=jp 254 : * 255 : * segdes xmatri 256 : * 257 : 10 continue 258 : if(nelri.ne.nelrig) then 259 : nelrig=nelri 260 : segadj xmatri 261 : endif 262 : * 263 : else 264 : * 265 : * creation des descriptifs communs a toutes les raideurs 266 : * 267 : imo=1 268 : if(ifour.eq.0) imo=3 269 : nligrp=5 270 : nligrd=nligrp 271 : segini descr 272 : irigel(3,1)=descr 273 : lisinc(1)='LX' 274 : lisdua(1)='FLX' 275 : noelep(1)=1 276 : noeled(1)=1 277 : do 70 i=2,5,2 278 : lisinc(i)=modepl(imo) 279 : lisinc(i+1)=modepl(imo+1) 280 : lisdua(i)=moforc(imo) 281 : lisdua(i+1)=moforc(imo+1) 282 : noelep(i)=(i+2)/2 283 : noelep(i+1)=noelep(i) 284 : noeled(i)=noelep(i) 285 : noeled(i+1)=noelep(i) 286 : 70 continue 287 : * 288 : * creation du segment imatri 289 : * 290 : nelrig=50 291 : nelri=0 292 : segini xmatri 293 : irigel(4,1)=xmatri 294 : * 295 : * ce qu'on cree est une égalité 296 : * 297 : irigel(6,1)=0 298 : * 299 : * ce qu'on cree est symetrique 300 : * 301 : irigel(7,1)=0 302 : * 303 : * boucle sur le maillage support 304 : * 305 : do 50 iel=1,ipt1.num(/2) 306 : * 307 : * noeud maitre 308 : * 309 : ip1=ipt1.num(2,iel) 310 : xp1=xcoor((ip1-1)*(idim+1)+1) 311 : yp1=xcoor((ip1-1)*(idim+1)+2) 312 : * 313 : * deuxieme noeud maitre fictif 314 : * 315 : if(iel.lt.ipt1.num(/2)) then 316 : ip2=ipt1.num(2,iel+1) 317 : jcm=1 318 : else 319 : ip2=ipt1.num(2,iel-1) 320 : jcm=-1 321 : endif 322 : xp2=xcoor((ip2-1)*(idim+1)+1) 323 : yp2=xcoor((ip2-1)*(idim+1)+2) 324 : xr=(xp2-xp1)*jcm 325 : yr=(yp2-yp1)*jcm 326 : * 327 : * noeud esclave 328 : * 329 : jp=ipt1.num(3,iel) 330 : xp=xcoor((jp-1)*(idim+1)+1) 331 : yp=xcoor((jp-1)*(idim+1)+2) 332 : * 333 : * verification que pas relation du point sur lui meme 334 : * 335 : if (jp.eq.ip1) goto 50 336 : * 337 : * detection du contact : point materiel 338 : * 339 : if (iret.eq.0) then 340 : if(((xp-xp1)*yr-(yp-yp1)*xr).gt.0.d0) goto 50 341 : vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2)) 342 : vn1=vec(1)/vnorm 343 : vn2=vec(2)/vnorm 344 : endif 345 : * 346 : * detection du contact : projectile a nez plat 347 : * 348 : if (iret.eq.1) then 349 : tv=vec(1)/vec(2) 350 : aa=(xp-xp1)*tv 351 : if(-(yp+aa-yp1)*xr.gt.0.d0) goto 50 352 : vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2)) 353 : if(abs(xp-xp1).gt.abs(vlarg*vec(2)/vnorm)) goto 50 354 : vn1=vec(1)/vnorm 355 : vn2=vec(2)/vnorm 356 : endif 357 : * 358 : * detection du contact : projectile a nez hemispherique 359 : * 360 : if (iret.eq.2) then 361 : vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2)) 362 : if(abs(xp-xp1).gt.vrayo*abs(vec(2))/vnorm) goto 50 363 : ypo=(vrayo/sqrt(1+(vec(1)*vec(1)/vec(2)*vec(2))))+yp 364 : xpo=(ypo-yp)*vec(1)/vec(2)+xp 365 : yy=ypo-sqrt(vrayo*vrayo-(xpo-xp1)*(xpo-xp1)) 366 : if(-(yy-yp1)*xr.gt.0.d0) goto 50 367 : vnn=sqrt((xpo-xp1)*(xpo-xp1)+(ypo-yp1)*(ypo-yp1)) 368 : vn1=(xp1-xpo)/vnn 369 : vn2=(yp1-ypo)/vnn 370 : endif 371 : * 372 : * detection du contact : projectile a nez conique 373 : * 374 : if (iret.eq.3) then 375 : ial=nint(sign(1.D0,(xp-xp1))) 376 : vangl2=vangl-ial*atan(vec(1)/vec(2)) 377 : aa=abs(xp-xp1)*tan(vangl2) 378 : if(-(yp+aa-yp1)*xr.gt.0.d0) goto 50 379 : vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2)) 380 : if(abs(xp-xp1).gt.vlarg*abs(vec(2))/vnorm) goto 50 381 : v1=cos(vangl)*vec(1)+ial*sin(vangl)*vec(2) 382 : v2=-ial*sin(vangl)*vec(1)+cos(vangl)*vec(2) 383 : vnn=sqrt(v1*v1+v2*v2) 384 : vn1=v1/vnn 385 : vn2=v2/vnn 386 : endif 387 : * 388 : * on cree le xmatri C 389 : * 390 : * segini xmatri 391 : * imatri=irigel(4,1) 392 : nelri=nelri+1 393 : if(nelri.gt.nelrig) then 394 : nelrig=nelrig+50 395 : segadj xmatri 396 : endif 397 : * imattt(nelrig)=xmatri 398 : * 399 : * on remplit le xmatri 400 : * 401 : * multiplicateur 402 : re(1,1,nelri)= 0.d0 403 : * 404 : * noeud maitre ip1 405 : re(2,1,nelri)= vn1 406 : re(3,1,nelri)=-vn2 407 : * 408 : * noeud esclave jp 409 : re(4,1,nelri)=-vn1 410 : re(5,1,nelri)= vn2 411 : * 412 : * test de conditionnement 413 : do 95 il=1,5 414 : if(abs(re(il,1,nelri)).lt.1D-6) re(il,1,nelri)=0.D0 415 : if(re(il,1,nelri).gt.(1.D0-1D-6)) re(il,1,nelri)=1.D0 416 : if(re(il,1,nelri).lt.(-1.D0+1D-6)) re(il,1,nelri)=-1.D0 417 : 95 continue 418 : * 419 : * transposition du xmatri C 420 : do 40 il=1,1 421 : do 401 ic=1,5 422 : re(il,ic,nelri)=re(ic,il,nelri) 423 : 401 continue 424 : 40 continue 425 : * 426 : * modification du meleme 427 : * 428 : nbelem=num(/2)+1 429 : nbnn=3 430 : nbsous=0 431 : nbref=0 432 : segadj meleme 433 : num(1,nbelem)=ipt1.num(1,iel) 434 : num(2,nbelem)=ip1 435 : num(3,nbelem)=jp 436 : * 437 : * segdes xmatri 438 : * 439 : 50 continue 440 : if(nelrig.ne.nelri) then 441 : nelrig=nelri 442 : segadj xmatri 443 : endif 444 : * 445 : endif 446 : * 447 : * desactivation des segments 448 : * 449 : segdes xmatri,ipt1,mrigid,meleme,descr 450 : * 451 : * renvoi du resultat 452 : * 453 : if(nbelem.eq.0) then 454 : irigel(6,1)=0 454 : call ecrobj('ANNULE',mrigid) 455 : else 455 : call ecrobj('RIGIDITE',mrigid) 456 : endif 457 : * 458 : end 459 : 460 :
© Cast3M 2003 - All rights reserved.
Disclaimer