C RELTUY SOURCE PV090527 26/04/30 21:16:11 12529 subroutine reltuy implicit real*8(a-h,o-z) implicit integer (i-n) -INC PPARAM -INC CCOPTIO -INC SMRIGID -INC SMELEME -INC SMCOORD -INC CCREEL segment icpr(nbpts) SEGMENT ISEG3 INTEGER NIZO(NxyZo+1) ENDSEGMENT * SEGMENT ISEG4 * INTEGER NUMZO(NxyZo) * ENDSEGMENT SEGMENT ISEG5 INTEGER NNMEL(ILON),IDEJ(nel),ipos(na+1) ENDSEGMENT SEGMENT ISEG6 INTEGER numpo(illo) ENDSEGMENT segment izon(na) DATA EPSI/1.D-4/ * * lecture des donnees on met le maillage de seg2 ou seg3 dans ipt1 * call lirobj( 'MAILLAGE',ipt1,1,iretou) if(ierr.ne.0) return call lirobj('MAILLAGE' ,ipt2,1,iretou) if(ierr.ne.0) return segact ipt1,ipt2 if(ipt1.lisous(/1).ne.0) then if(ipt2.itypel.ne.2.and.ipt2.itypel.ne.3)then call erreur (19) return endif ipta=ipt2 ipt2=ipt1 ipt1=ipta else if (ipt1.itypel.ne.2.and.ipt1.itypel.ne.3) then if(ipt2.itypel.ne.2.and.ipt2.itypel.ne.3)then call erreur(19) return else ipta=ipt2 ipt2=ipt1 ipt1=ipta endif endif endif call change (ipt2,1) if( ipt1.itypel.eq.3) call change (ipt1,2) if( ipt1.lisous(/1).ne.0) then * on suppose un seul type d'elements pour la tuyauterir call erreur(19) return endif npipt2=ipt2.num(/2) * * * lecture facultative de la plus grande section d'un tuyau * call lirree(sect,1,iretou) if(ierr.ne.0) return * *pour zoner on prend la longueurn max d'un element de ajouter a sect * xlon=0.d0 ider=ipt1.num(/1) do iel=1,ipt1.num(/2) ia= ipt1.num(1,iel) ib=ipt1.num(ider,iel) ia1=ia*(idim+1)-idim xa=xcoor(ia1) ya=xcoor ( ia1+1) za=xcoor ( ia1+2) ib1=ib*(idim+1)-idim xb=xcoor(ib1) yb=xcoor ( ib1+1) zb=xcoor ( ib1+2) xlo= sqrt( (xa-xb)*(xa-xb)+(ya-yb)*(ya-yb)+(za-zb)*(za-zb)) if( xlo . gt. xlon) xlon=xlo enddo xpre=epsi*xlon xcri= sqrt(xlon*xlon + sect*sect) * * * zonage on cherche les min et max suivant les trois directions et la taille *max des elements de la ligne de tuyauterie suivant les trois directions * * ITR = 0 XZO=0.D0 YZO=0.D0 ZZO=0.D0 XZA=XGRAND YZA=XGRAND ZZA=XGRAND XMI=XGRAND XMA=-XGRAND YMI=XGRAND YMA=-XGRAND if(idim.eq.3) then ZMI=XGRAND ZMA=-XGRAND else zmi=0.d0 zma=0.d0 endif meleme=ipt1 ipasag=0 idim1=idim+1 254 ipasag=ipasag+1 if(ipasag.eq.1) then meleme=ipt1 else meleme=ipt2 endif nbnn= num(/1) nbel=num(/2) do 200 i1=1,nbel DO 201 I2 = 1,NBNN IB=NUM(I2,I1) IA=(IB-1)*IDIM1 IF(XCOOR(IA+1).LT.XMI) XMI=XCOOR(IA+1) IF(XCOOR(IA+1).GT.XMA) XMA=XCOOR(IA+1) IF(XCOOR(IA+2).LT.YMI) YMI=XCOOR(IA+2) IF(XCOOR(IA+2).GT.YMA) YMA=XCOOR(IA+2) IF ( IDIM.EQ.3 ) THEN IF(XCOOR(IA+3).LT.ZMI) ZMI=XCOOR(IA+3) IF(XCOOR(IA+3).GT.ZMA) ZMA=XCOOR(IA+3) ENDIF 201 CONTINUE 200 CONTINUE If(ipasag.eq.1) go to 254 nzo=1 xmi= xmi-( xma-xmi ) / 1.d6 ymi= ymi-( yma-ymi ) / 1.d6 xma=xma + ( xma-xmi ) / 1.d6 yma=yma + ( yma-ymi ) / 1.d6 nxo= int(( xma - xmi ) / xcri) +1 nyo= int(( yma - ymi ) / xcri) + 1 if(idim.eq.3) then zmi= zmi-( zma-zmi ) / 1.d6 zma=zma+( zma-zmi ) / 1.d6 nzo= int (( zma - zmi ) / xcri) + 1 endif nxyo=nxo*nyo nxyzo= nxo*nyo*nzo * write(6,*) ' xcri nxo nyo nzo nxyo nxyzo' * write(6,*) xcri, nxo ,nyo, nzo, nxyo ,nxyzo * write(6,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma * * on fait une numerotation locale pour les points de la ligne * on ne s'occupe que des points extremes * na=0 segini icpr do iel=1,ipt1.num(/2) do ip=1,ipt1.num(/1) ia= ipt1.num(ip,iel) if(icpr(ia).eq.0) then na=na+1 icpr(ia)=na endif enddo enddo * write(6,*) ' na nxyo nxo' , na,nxyo , nxo * * on cherche dans quelle zone est chacun des points de la tuyauterie * et on dresse le tableau la liste des points par zone segini iseg3,izon zp=0.d0 do iu=1,icpr(/1) if(icpr(iu).ne.0) then ip=icpr(iu) xp= xcoor( (iu-1)*(idim+1) +1) yp=xcoor( (iu-1)*(idim+1) +2) if( idim.eq.3)ZP= xcoor((iu-1)*(idim+1) +3) nn= (int(( zp-zmi)/xcri) )*nxyo $ +(int(( yp-ymi)/xcri) )*nxo $ +(int(( xp-xmi)/xcri) + 1) nizo(nn)=nizo(nn)+1 izon(ip)=nn * write(6,*) 'xp yp zp iu ip izon(ip) ',xp,yp,zp ,iu,ip,izon(ip) endif enddo do ia=2,nxyzo+1 nizo(ia)=nizo(ia)+nizo(ia-1) enddo * write(6,*) ' nizo ' , (nizo(iou),iou=1,nizo(/1)) illo=nizo(nxyzo+1) segini iseg6 do ip1=1,icpr(/1) inp=icpr(ip1) if( inp.ne.0) then nn=izon(inp) ipla=nizo(nn) nizo(nn)=nizo(nn)-1 numpo(ipla)=ip1 endif enddo * write(6,*) ' nizo ' , (nizo(iou),iou=1,nizo(/1)) * write(6,*) ' numpo ',(numpo(iou),iou=1,numpo(/1)) * on construit le tableua donnant le numero des elements touchant un noeud * en numerotation locale ilon=2*na nel=ipt1.num(/2) segini iseg5 do ia=1,ipt1.num(/2) do ib=1,ipt1.num(/1) ip1=ipt1.num(ib,ia) iy=icpr(ip1) ipos(iy)=ipos(iy)+1 enddo enddo do ia=2,na+1 ipos(ia)=ipos(ia-1)+ipos(ia) enddo do ia=1,ipt1.num(/2) do ib=1,ipt1.num(/1) ip1=ipt1.num(ib,ia) iy=icpr(ip1) ipla=ipos(iy) ipos(iy)=ipos(iy)-1 nnmel(ipla)=ia enddo enddo * write(6,*) ' ipos ' , ( ipos(iou),iou=1,ipos(/1)) * write(6,*) ' nnmel ' , ( nnmel(iou),iou=1,nnmel(/1)) * * on cree tout de suite les points des multiplicateur * segact mcoord*mod iposmu=nbpts nbpts=iposmu + npipt2 segadj mcoord * * pour les relations il peut y avoir des relations sur tout le segment ou sur * seulement un point. ON cree les deux types, on ajustera plus tard. * nbnn= ipt1.num(/1)+2 nbelem= ipt2.num(/2) nbsous=0 nbref=0 segini ipt4 ipt4.itypel=22 nrigel=2 segini mrigid mtymat='RIGIDITE' iforig=ifour coerig(1)=1.d0 coerig(2)=1.d0 irigel(1,1)=ipt4 nligrp=nbnn nligrd=nbnn segini descr irigel(3,1)=descr lisinc(1)='LX' lisdua(1)='FLX' noelep(1)=1 noeleD(1)=1 do iu=2,nbnn noelep(iu)=iu noeled(iu)=iu lisinc(iu)='T' lisdua(iu)='Q' enddo segdes descr nelrig= nbelem rigrel=0 segini xmatr4 irigel(4,1)= xmatr4 irigel(5,1)=nifour nel4=0 nel5=0 nbnn=3 * write(6,*) ' à la creation de ipt5 nbelem ' , nbelem segini ipt5 ipt5.itypel=22 irigel(1,2)=ipt5 nligrp=nbnn nligrd=nbnn rigrel=0 segini xmatr5 irigel(4,2)=xmatr5 irigel(5,2)=nifour segini descr irigel(3,2)=descr lisinc(1)='LX' lisdua(1)='FLX' noelep(1)=1 noeleD(1)=1 noelep(2)=2 noeled(2)=2 noelep(3)=3 noeled(3)=3 lisinc(2)='T' lisinc(3)='T' lisdua(2)='Q' lisdua(3)='Q' segdes descr nbnn=ipt1.num(/1) nz=1 do ipp=1,ipt2.num(/2) ip= ipt2.num(1,ipp) iqv=0 xp= xcoor( (ip-1)*(idim+1) +1) yp=xcoor( (ip-1)*(idim+1) +2) zp=0.d0 if(idim.eq.3)ZP= xcoor((ip-1)*(idim+1) +3) nx= int(( xp-xmi)/xcri) + 1 ny= int(( yp-ymi)/xcri) + 1 nz= int(( zp-zmi)/xcri) + 1 nn= (nz-1)*nxyo +(ny-1)*nxo + nx * write(6,*) ' point ' , ip , ' zonze ' , nn, xp, yp,zp,nzo * write(6,*) ' point , ip ,nx ny nz ' , nx , ny , nz ield=0 xdi=xgrand nnxmi= max ( nx-1,1) nnxma= min (nx+1,nxo) nnymi= max(1,ny-1) nnyma= min(nyo,nY+1) nnzmi= max(1,nz-1) nnzma= min(nzo,nz+1) * write(6,*) 'mima', nnxmi,nnxma,nnymi,nnyma,nnzmi,nnzma do ix=nnxmi,nnxma do iy=nnymi,nnyma do iz=nnzmi,nnzma nnes=(iz-1)*nxyo + (iy-1)*nxo+ix itp =nizo(nnes+1)-nizo(nnes) * write(6,* ) ' ip ix iy iz nnes itp' , ix,iy,iz, nnes,itp if (itp.ne.0) then idec= nizo(nnes) do im=1,itp ippp=numpo( idec + im) ipl=icpr(ippp) * write(6,*) ' ippp ipl ' , ippp,ipl iplael=ipos(ipl)+1 do ielt=iplael,ipos(ipl+1) iel= nnmel(ielt) * write(6,*) ' ielt iel nbnn' , ielt,iel,nbnn ipa=(ipt1.num(1,iel)-1)*(idim+1) ipb=(ipt1.num(nbnn,iel)-1)*(idim+1) * write(6,*) ' ipa ipb ' , ipa/(idim+1), ipb/(idim+1) xa = xcoor(ipa+1) ya = xcoor(ipa+2) xb = xcoor(ipb+1) yb = xcoor(ipb+2) if(idim.eq.3) then za = xcoor(ipa+3) zb = xcoor(ipb+3) else za=0.d0 zb=0.d0 endif * * calcul de la distance de ip a la droit A B et verif si pt interne ou non * xab=xb-xa yab=yb-ya zab=zb-za xlo= sqrt(xab*xab + yab * yab + zab * zab) xab=xab/xlo yab=yab/xlo zab=zab/xlo xloi= xab*( xp-xa)+yab*(yp-ya)+zab*(zp-za) xloi=xloi/xlo * if( ip.eq.118476) then * write(6,*) ' xloi xpre xlo ' , xloi, xpre,xlo * write(6,*)' ipa ipb ',ipt1.num(1,iel),ipt1.num(nbnn,iel) * write(6,*) ' xab yab zab ',xab,yab,zab * write(6,*) ' xa ya za xp yp zp',xa,ya,za,xp,yp,zp * endif if( xloi.ge.epsi.and. xloi.le.(1.D0-epsi) ) then vx= yab*(zp-za) - zab*(yp-ya) vy= zab*(xp-xa) - xab*(zp-za) vz= xab*(yp-ya) - yab*(xp-xa) xdis=sqrt(vx*vx + vy*vy + vz*vz) * if( ip.eq.118476) then * write(6,*) ' on a trouve un candidat xdis' , xdis * endif if( xdis.lt.xdi) then ield=iel xdi=xdis xrap=xloi endif endif enddo enddo endif enddo enddo enddo * a-t-on trouvé un candidat? if(ield.ne.0) then iposmu=iposmu+1 nel4=nel4+1 ipt4.num(1,nel4)=iposmu ipt4.num(2,nel4)=ip ipt4.num(3,nel4)=ipt1.num(1,ield) ipt4.num(4,nel4)=ipt1.num(2,ield) xmatr4.re(1,2,nel4)=1.d0 xmatr4.re(1,3,nel4)= xrap -1.D0 xmatr4.re(1,4,nel4)= -xrap xmatr4.re(2,1,nel4)=1.d0 xmatr4.re(3,1,nel4)= xrap -1.D0 xmatr4.re(4,1,nel4)= -xrap else * il faut recommencer pour trouver le point le plus proche et ecrire une * relation sur un seul point xdi=xgrand * write(6,*) 'nnxmi,nnxma',nnxmi,nnxma * write(6,*) 'nnymi,nnyma',nnymi,nnyma * write(6,*) 'nnzmi,nnzma',nnzmi,nnzma do ix=nnxmi,nnxma do iy=nnymi,nnyma do iz=nnzmi,nnzma nnes= (iz-1)*nxyo + (iy-1)*nxo+ix * write(6,*) 'nizo(nnes)+1,nizo(nnes+1)',nizo(nnes),nizo(nnes+1) do iqz=nizo(nnes)+1,nizo(nnes+1) iq=numpo(iqz) iqq=(iq-1)*(idim+1) xq=xcoor(iqq+1) yq=xcoor(iqq+2) zq=0.d0 if( idim.eq.3) zq=xcoor(iqq+3) xdis=(xp-xq)*(xp-xq) + (yp-yq)*(yp-yq) + ( zp-zq)*(zp-zq) * write(6,*) ' point essaye iq xdis ' , iq ,xdis , xdi if( xdi.gt.xdis) then xdi=xdis iqv=iq endif enddo enddo enddo enddo nel5=nel5+1 iposmu=iposmu+1 ipt5.num(1,nel5)=iposmu ipt5.num(2,nel5)=ip ipt5.num(3,nel5)=iqv xmatr5.re(1,2,nel5)=1.d0 xmatr5.re(2,1,nel5)=1.D0 xmatr5.re(1,3,nel5)=-1.d0 xmatr5.re(3,1,nel5)=-1.d0 endif enddo * * ajustement des objets rigidité * * write(6,*) ' nbelem nel4 nel5 ' , nbelem,nel4 , nel5 * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8) if(nel4.ne.0) then if(ipt4.num(/2).ne.nel4) then nelrig=nel4 nbelem=nel4 nbnn= ipt4.num(/1) segadjipt4 nligrp=xmatr4.re(/1) nligrd=xmatr4.re(/2) rigrel=0 segadj xmatr4 nelrig= nel5 nbelem=nel5 nbnn=ipt5.num(/1) segadj ipt5 nligrp=xmatr5.re(/1) nligrd=xmatr5.re(/2) rigrel=0 segadj xmatr5 else * write(6,*) ' on passe bien la' nrigel=1 segadj mrigid * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8) endif else do io=1,irigel(/1) irigel(io,1)=irigel(io,2) enddo nrigel=1 segadj mrigid endif * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8) segdes ipt4,ipt5,xmatr4,xmatr5,mrigid segsup izon,iseg3,iseg5,iseg6,icpr call ecrobj('RIGIDITE',mrigid) return end