1 : C RELTUY SOURCE CHAT 12/09/11 21:15:11 7499 2 : subroutine reltuy 3 : implicit real*8(a-h,o-z) 4 : implicit integer (i-n) 5 : -INC CCOPTIO 6 : -INC SMRIGID 7 : -INC SMELEME 8 : -INC SMCOORD 9 : -INC CCREEL 10 : segment icpr(xcoor(/1)/(idim+1)) 11 : SEGMENT ISEG3 12 : INTEGER NIZO(NxyZo+1) 13 : ENDSEGMENT 14 : * SEGMENT ISEG4 15 : * INTEGER NUMZO(NxyZo) 16 : * ENDSEGMENT 17 : SEGMENT ISEG5 18 : INTEGER NNMEL(ILON),IDEJ(nel),ipos(na+1) 19 : ENDSEGMENT 20 : SEGMENT ISEG6 21 : INTEGER numpo(illo) 22 : ENDSEGMENT 23 : segment izon(na) 24 : DATA EPSI/1.D-4/ 25 : 26 : * 27 : * lecture des donnees on met le maillage de seg2 ou seg3 dans ipt1 28 : * 28 : call lirobj( 'MAILLAGE',ipt1,1,iretou) 29 : if(ierr.ne.0) return 29 : call lirobj('MAILLAGE' ,ipt2,1,iretou) 30 : if(ierr.ne.0) return 31 : segact ipt1,ipt2 32 : if(ipt1.lisous(/1).ne.0) then 33 : if(ipt2.itypel.ne.2.and.ipt2.itypel.ne.3)then 33 : call erreur (19) 34 : return 35 : endif 36 : ipta=ipt2 37 : ipt2=ipt1 38 : ipt1=ipta 39 : else 40 : if (ipt1.itypel.ne.2.and.ipt1.itypel.ne.3) then 41 : if(ipt2.itypel.ne.2.and.ipt2.itypel.ne.3)then 41 : call erreur(19) 42 : return 43 : else 44 : ipta=ipt2 45 : ipt2=ipt1 46 : ipt1=ipta 47 : endif 48 : endif 49 : endif 49 : call change (ipt2,1) 49 : if( ipt1.itypel.eq.3) call change (ipt1,2) 50 : if( ipt1.lisous(/1).ne.0) then 51 : * on suppose un seul type d'elements pour la tuyauterir 51 : call erreur(19) 52 : return 53 : endif 54 : npipt2=ipt2.num(/2) 55 : * 56 : * 57 : * lecture facultative de la plus grande section d'un tuyau 58 : * 58 : call lirree(sect,1,iretou) 59 : if(ierr.ne.0) return 60 : * 61 : *pour zoner on prend la longueurn max d'un element de ajouter a sect 62 : * 63 : xlon=0.d0 64 : ider=ipt1.num(/1) 65 : do iel=1,ipt1.num(/2) 66 : ia= ipt1.num(1,iel) 67 : ib=ipt1.num(ider,iel) 68 : ia1=ia*(idim+1)-idim 69 : xa=xcoor(ia1) 70 : ya=xcoor ( ia1+1) 71 : za=xcoor ( ia1+2) 72 : ib1=ib*(idim+1)-idim 73 : xb=xcoor(ib1) 74 : yb=xcoor ( ib1+1) 75 : zb=xcoor ( ib1+2) 76 : xlo= sqrt( (xa-xb)*(xa-xb)+(ya-yb)*(ya-yb)+(za-zb)*(za-zb)) 77 : if( xlo . gt. xlon) xlon=xlo 78 : enddo 79 : xpre=epsi*xlon 80 : xcri= sqrt(xlon*xlon + sect*sect) 81 : * 82 : * 83 : * zonage on cherche les min et max suivant les trois directions et la taille 84 : *max des elements de la ligne de tuyauterie suivant les trois directions 85 : * 86 : * 87 : ITR = 0 88 : XZO=0.D0 89 : YZO=0.D0 90 : ZZO=0.D0 91 : XZA=XGRAND 92 : YZA=XGRAND 93 : ZZA=XGRAND 94 : XMI=XGRAND 95 : XMA=-XGRAND 96 : YMI=XGRAND 97 : YMA=-XGRAND 98 : if(idim.eq.3) then 99 : ZMI=XGRAND 100 : ZMA=-XGRAND 101 : else 102 : zmi=0.d0 103 : zma=0.d0 104 : endif 105 : meleme=ipt1 106 : ipasag=0 107 : idim1=idim+1 108 : 254 ipasag=ipasag+1 109 : if(ipasag.eq.1) then 110 : meleme=ipt1 111 : else 112 : meleme=ipt2 113 : endif 114 : nbnn= num(/1) 115 : nbel=num(/2) 116 : do 200 i1=1,nbel 117 : DO 201 I2 = 1,NBNN 118 : IB=NUM(I2,I1) 119 : IA=(IB-1)*IDIM1 120 : IF(XCOOR(IA+1).LT.XMI) XMI=XCOOR(IA+1) 121 : IF(XCOOR(IA+1).GT.XMA) XMA=XCOOR(IA+1) 122 : IF(XCOOR(IA+2).LT.YMI) YMI=XCOOR(IA+2) 123 : IF(XCOOR(IA+2).GT.YMA) YMA=XCOOR(IA+2) 124 : IF ( IDIM.EQ.3 ) THEN 125 : IF(XCOOR(IA+3).LT.ZMI) ZMI=XCOOR(IA+3) 126 : IF(XCOOR(IA+3).GT.ZMA) ZMA=XCOOR(IA+3) 127 : ENDIF 128 : 201 CONTINUE 129 : 200 CONTINUE 130 : If(ipasag.eq.1) go to 254 131 : nzo=1 132 : xmi= xmi-( xma-xmi ) / 1.d6 133 : ymi= ymi-( yma-ymi ) / 1.d6 134 : xma=xma + ( xma-xmi ) / 1.d6 135 : yma=yma + ( yma-ymi ) / 1.d6 136 : nxo= int(( xma - xmi ) / xcri) +1 137 : nyo= int(( yma - ymi ) / xcri) + 1 138 : if(idim.eq.3) then 139 : zmi= zmi-( zma-zmi ) / 1.d6 140 : zma=zma+( zma-zmi ) / 1.d6 141 : nzo= int (( zma - zmi ) / xcri) + 1 142 : endif 143 : nxyo=nxo*nyo 144 : nxyzo= nxo*nyo*nzo 145 : * write(6,*) ' xcri nxo nyo nzo nxyo nxyzo' 146 : * write(6,*) xcri, nxo ,nyo, nzo, nxyo ,nxyzo 147 : * write(6,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma 148 : * 149 : * on fait une numerotation locale pour les points de la ligne 150 : * on ne s'occupe que des points extremes 151 : * 152 : na=0 153 : segini icpr 154 : do iel=1,ipt1.num(/2) 155 : do ip=1,ipt1.num(/1) 156 : ia= ipt1.num(ip,iel) 157 : if(icpr(ia).eq.0) then 158 : na=na+1 159 : icpr(ia)=na 160 : endif 161 : enddo 162 : enddo 163 : * write(6,*) ' na nxyo nxo' , na,nxyo , nxo 164 : * 165 : * on cherche dans quelle zone est chacun des points de la tuyauterie 166 : * et on dresse le tableau la liste des points par zone 167 : 168 : segini iseg3,izon 169 : zp=0.d0 170 : do iu=1,icpr(/1) 171 : if(icpr(iu).ne.0) then 172 : ip=icpr(iu) 173 : xp= xcoor( (iu-1)*(idim+1) +1) 174 : yp=xcoor( (iu-1)*(idim+1) +2) 175 : if( idim.eq.3)ZP= xcoor((iu-1)*(idim+1) +3) 176 : nn= (int(( zp-zmi)/xcri) )*nxyo 177 : $ +(int(( yp-ymi)/xcri) )*nxo 178 : $ +(int(( xp-xmi)/xcri) + 1) 179 : nizo(nn)=nizo(nn)+1 180 : izon(ip)=nn 181 : * write(6,*) 'xp yp zp iu ip izon(ip) ',xp,yp,zp ,iu,ip,izon(ip) 182 : endif 183 : enddo 184 : do ia=2,nxyzo+1 185 : nizo(ia)=nizo(ia)+nizo(ia-1) 186 : enddo 187 : * write(6,*) ' nizo ' , (nizo(iou),iou=1,nizo(/1)) 188 : illo=nizo(nxyzo+1) 189 : segini iseg6 190 : do ip1=1,icpr(/1) 191 : inp=icpr(ip1) 192 : if( inp.ne.0) then 193 : nn=izon(inp) 194 : ipla=nizo(nn) 195 : nizo(nn)=nizo(nn)-1 196 : numpo(ipla)=ip1 197 : endif 198 : enddo 199 : * write(6,*) ' nizo ' , (nizo(iou),iou=1,nizo(/1)) 200 : * write(6,*) ' numpo ',(numpo(iou),iou=1,numpo(/1)) 201 : * on construit le tableua donnant le numero des elements touchant un noeud 202 : * en numerotation locale 203 : ilon=2*na 204 : nel=ipt1.num(/2) 205 : segini iseg5 206 : do ia=1,ipt1.num(/2) 207 : do ib=1,ipt1.num(/1) 208 : ip1=ipt1.num(ib,ia) 209 : iy=icpr(ip1) 210 : ipos(iy)=ipos(iy)+1 211 : enddo 212 : enddo 213 : do ia=2,na+1 214 : ipos(ia)=ipos(ia-1)+ipos(ia) 215 : enddo 216 : do ia=1,ipt1.num(/2) 217 : do ib=1,ipt1.num(/1) 218 : ip1=ipt1.num(ib,ia) 219 : iy=icpr(ip1) 220 : ipla=ipos(iy) 221 : ipos(iy)=ipos(iy)-1 222 : nnmel(ipla)=ia 223 : enddo 224 : enddo 225 : * write(6,*) ' ipos ' , ( ipos(iou),iou=1,ipos(/1)) 226 : * write(6,*) ' nnmel ' , ( nnmel(iou),iou=1,nnmel(/1)) 227 : * 228 : * on cree tout de suite les points des multiplicateur 229 : * 230 : iposmu=xcoor(/1)/(idim+1) 231 : nbpts=iposmu + npipt2 232 : segadj mcoord 233 : 234 : * 235 : * pour les relations il peut y avoir des relations sur tout le segment ou sur 236 : * seulement un point. ON cree les deux types, on ajustera plus tard. 237 : * 238 : nbnn= ipt1.num(/1)+2 239 : nbelem= ipt2.num(/2) 240 : nbsous=0 241 : nbref=0 242 : segini ipt4 243 : ipt4.itypel=22 244 : nrigel=2 245 : segini mrigid 246 : mtymat='RIGIDITE' 247 : coerig(1)=1.d0 248 : coerig(2)=1.d0 249 : irigel(1,1)=ipt4 250 : nligrp=nbnn 251 : nligrd=nbnn 252 : segini descr 253 : irigel(3,1)=descr 254 : lisinc(1)='LX' 255 : lisdua(1)='FLX' 256 : noelep(1)=1 257 : noeleD(1)=1 258 : do iu=2,nbnn 259 : noelep(iu)=iu 260 : noeled(iu)=iu 261 : lisinc(iu)='T' 262 : lisdua(iu)='Q' 263 : enddo 264 : segdes descr 265 : nelrig= nbelem 266 : segini xmatr4 267 : irigel(4,1)= xmatr4 268 : irigel(5,1)=nifour 269 : nel4=0 270 : nel5=0 271 : nbnn=3 272 : * write(6,*) ' à la creation de ipt5 nbelem ' , nbelem 273 : segini ipt5 274 : ipt5.itypel=22 275 : irigel(1,2)=ipt5 276 : nligrp=nbnn 277 : nligrd=nbnn 278 : segini xmatr5 279 : irigel(4,2)=xmatr5 280 : irigel(5,2)=nifour 281 : segini descr 282 : irigel(3,2)=descr 283 : lisinc(1)='LX' 284 : lisdua(1)='FLX' 285 : noelep(1)=1 286 : noeleD(1)=1 287 : noelep(2)=2 288 : noeled(2)=2 289 : noelep(3)=3 290 : noeled(3)=3 291 : lisinc(2)='T' 292 : lisinc(3)='T' 293 : lisdua(2)='Q' 294 : lisdua(3)='Q' 295 : segdes descr 296 : nbnn=ipt1.num(/1) 297 : nz=1 298 : do ipp=1,ipt2.num(/2) 299 : ip= ipt2.num(1,ipp) 300 : iqv=0 301 : xp= xcoor( (ip-1)*(idim+1) +1) 302 : yp=xcoor( (ip-1)*(idim+1) +2) 303 : zp=0.d0 304 : if(idim.eq.3)ZP= xcoor((ip-1)*(idim+1) +3) 305 : nx= int(( xp-xmi)/xcri) + 1 306 : ny= int(( yp-ymi)/xcri) + 1 307 : nz= int(( zp-zmi)/xcri) + 1 308 : nn= (nz-1)*nxyo +(ny-1)*nxo + nx 309 : * write(6,*) ' point ' , ip , ' zonze ' , nn, xp, yp,zp,nzo 310 : * write(6,*) ' point , ip ,nx ny nz ' , nx , ny , nz 311 : ield=0 312 : xdi=xgrand 313 : nnxmi= max ( nx-1,1) 314 : nnxma= min (nx+1,nxo) 315 : nnymi= max(1,ny-1) 316 : nnyma= min(nyo,nY+1) 317 : nnzmi= max(1,nz-1) 318 : nnzma= min(nzo,nz+1) 319 : * write(6,*) 'mima', nnxmi,nnxma,nnymi,nnyma,nnzmi,nnzma 320 : do ix=nnxmi,nnxma 321 : do iy=nnymi,nnyma 322 : do iz=nnzmi,nnzma 323 : nnes=(iz-1)*nxyo + (iy-1)*nxo+ix 324 : itp =nizo(nnes+1)-nizo(nnes) 325 : * write(6,* ) ' ip ix iy iz nnes itp' , ix,iy,iz, nnes,itp 326 : if (itp.ne.0) then 327 : idec= nizo(nnes) 328 : do im=1,itp 329 : ippp=numpo( idec + im) 330 : ipl=icpr(ippp) 331 : * write(6,*) ' ippp ipl ' , ippp,ipl 332 : iplael=ipos(ipl)+1 333 : do ielt=iplael,ipos(ipl+1) 334 : iel= nnmel(ielt) 335 : * write(6,*) ' ielt iel nbnn' , ielt,iel,nbnn 336 : ipa=(ipt1.num(1,iel)-1)*(idim+1) 337 : ipb=(ipt1.num(nbnn,iel)-1)*(idim+1) 338 : * write(6,*) ' ipa ipb ' , ipa/(idim+1), ipb/(idim+1) 339 : xa = xcoor(ipa+1) 340 : ya = xcoor(ipa+2) 341 : xb = xcoor(ipb+1) 342 : yb = xcoor(ipb+2) 343 : if(idim.eq.3) then 344 : za = xcoor(ipa+3) 345 : zb = xcoor(ipb+3) 346 : else 347 : za=0.d0 348 : zb=0.d0 349 : endif 350 : * 351 : * calcul de la distance de ip a la droit A B et verif si pt interne ou non 352 : * 353 : xab=xb-xa 354 : yab=yb-ya 355 : zab=zb-za 356 : xlo= sqrt(xab*xab + yab * yab + zab * zab) 357 : xab=xab/xlo 358 : yab=yab/xlo 359 : zab=zab/xlo 360 : xloi= xab*( xp-xa)+yab*(yp-ya)+zab*(zp-za) 361 : xloi=xloi/xlo 362 : * if( ip.eq.118476) then 363 : * write(6,*) ' xloi xpre xlo ' , xloi, xpre,xlo 364 : * write(6,*)' ipa ipb ',ipt1.num(1,iel),ipt1.num(nbnn,iel) 365 : * write(6,*) ' xab yab zab ',xab,yab,zab 366 : * write(6,*) ' xa ya za xp yp zp',xa,ya,za,xp,yp,zp 367 : * endif 368 : if( xloi.ge.epsi.and. xloi.le.(1.D0-epsi) ) then 369 : vx= yab*(zp-za) - zab*(yp-ya) 370 : vy= zab*(xp-xa) - xab*(zp-za) 371 : vz= xab*(yp-ya) - yab*(xp-xa) 372 : xdis=sqrt(vx*vx + vy*vy + vz*vz) 373 : * if( ip.eq.118476) then 374 : * write(6,*) ' on a trouve un candidat xdis' , xdis 375 : * endif 376 : if( xdis.lt.xdi) then 377 : ield=iel 378 : xdi=xdis 379 : xrap=xloi 380 : endif 381 : endif 382 : enddo 383 : enddo 384 : endif 385 : enddo 386 : enddo 387 : enddo 388 : * a-t-on trouvé un candidat? 389 : if(ield.ne.0) then 390 : iposmu=iposmu+1 391 : nel4=nel4+1 392 : ipt4.num(1,nel4)=iposmu 393 : ipt4.num(2,nel4)=ip 394 : ipt4.num(3,nel4)=ipt1.num(1,ield) 395 : ipt4.num(4,nel4)=ipt1.num(2,ield) 396 : xmatr4.re(1,2,nel4)=1.d0 397 : xmatr4.re(1,3,nel4)= xrap -1.D0 398 : xmatr4.re(1,4,nel4)= -xrap 399 : xmatr4.re(2,1,nel4)=1.d0 400 : xmatr4.re(3,1,nel4)= xrap -1.D0 401 : xmatr4.re(4,1,nel4)= -xrap 402 : else 403 : * il faut recommencer pour trouver le point le plus proche et ecrire une 404 : * relation sur un seul point 405 : xdi=xgrand 406 : * write(6,*) 'nnxmi,nnxma',nnxmi,nnxma 407 : * write(6,*) 'nnymi,nnyma',nnymi,nnyma 408 : * write(6,*) 'nnzmi,nnzma',nnzmi,nnzma 409 : do ix=nnxmi,nnxma 410 : do iy=nnymi,nnyma 411 : do iz=nnzmi,nnzma 412 : nnes= (iz-1)*nxyo + (iy-1)*nxo+ix 413 : * write(6,*) 'nizo(nnes)+1,nizo(nnes+1)',nizo(nnes),nizo(nnes+1) 414 : do iqz=nizo(nnes)+1,nizo(nnes+1) 415 : iq=numpo(iqz) 416 : 417 : iqq=(iq-1)*(idim+1) 418 : xq=xcoor(iqq+1) 419 : yq=xcoor(iqq+2) 420 : zq=0.d0 421 : if( idim.eq.3) zq=xcoor(iqq+3) 422 : xdis=(xp-xq)*(xp-xq) + (yp-yq)*(yp-yq) + ( zp-zq)*(zp-zq) 423 : * write(6,*) ' point essaye iq xdis ' , iq ,xdis , xdi 424 : if( xdi.gt.xdis) then 425 : xdi=xdis 426 : iqv=iq 427 : endif 428 : enddo 429 : enddo 430 : enddo 431 : enddo 432 : nel5=nel5+1 433 : iposmu=iposmu+1 434 : ipt5.num(1,nel5)=iposmu 435 : ipt5.num(2,nel5)=ip 436 : ipt5.num(3,nel5)=iqv 437 : xmatr5.re(1,2,nel5)=1.d0 438 : xmatr5.re(2,1,nel5)=1.D0 439 : xmatr5.re(1,3,nel5)=-1.d0 440 : xmatr5.re(3,1,nel5)=-1.d0 441 : endif 442 : enddo 443 : * 444 : * ajustement des objets rigidité 445 : * 446 : 447 : * write(6,*) ' nbelem nel4 nel5 ' , nbelem,nel4 , nel5 448 : * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8) 449 : if(nel4.ne.0) then 450 : if(ipt4.num(/2).ne.nel4) then 451 : nelrig=nel4 452 : nbelem=nel4 453 : nbnn= ipt4.num(/1) 454 : segadjipt4 455 : nligrp=xmatr4.re(/1) 456 : nligrd=xmatr4.re(/2) 457 : segadj xmatr4 458 : nelrig= nel5 459 : nbelem=nel5 460 : nbnn=ipt5.num(/1) 461 : segadj ipt5 462 : nligrp=xmatr5.re(/1) 463 : nligrd=xmatr5.re(/2) 464 : segadj xmatr5 465 : else 466 : * write(6,*) ' on passe bien la' 467 : nrigel=1 468 : segadj mrigid 469 : * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8) 470 : endif 471 : else 472 : do io=1,irigel(/1) 473 : irigel(io,1)=irigel(io,2) 474 : enddo 475 : nrigel=1 476 : segadj mrigid 477 : endif 478 : * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8) 479 : segdes ipt4,ipt5,xmatr4,xmatr5,mrigid 480 : segsup izon,iseg3,iseg5,iseg6,icpr 480 : call ecrobj('RIGIDITE',mrigid) 481 : return 482 : end 483 : 484 :
© Cast3M 2003 - All rights reserved.
Disclaimer