Télécharger la plaquette Cast3M


fr en

Télécharger reltuy.eso

Retour à la liste

Numérotation des lignes : 

reltuy

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales