Télécharger la plaquette Cast3M


fr en

Télécharger redcon.eso

Retour à la liste

Numérotation des lignes : 

redcon

  1. C REDCON    SOURCE    CB215821  24/04/12    21:17:05     11897          
  2.       subroutine redcon(mmodel,mmode2)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.  
  6. -INC PPARAM
  7. -INC CCOPTIO
  8. -INC SMCOORD
  9. -INC SMELEME
  10. -INC CCREEL
  11. -INC SMMODEL
  12.       character*4 modepl(3),moforc(3)
  13.       data modepl /'UX  ','UY  ','UZ  '/
  14.       data moforc /'FX  ','FY  ','FZ  '/
  15. *
  16. *  lecture  du maillage support des conditions
  17. *
  18.       segact mmodel
  19.       segini,mmode2=mmodel
  20.       imoin=0
  21.       do  1 imo=1,kmodel(/1)
  22.         imodel=kmodel(imo)
  23.         segact imodel
  24.         do ifo=1,formod(/2)
  25.           if( formod(ifo).eq.'CONTACT' ) go to 2
  26.         enddo
  27.         mmode2.kmodel(imo-imoin)=imodel
  28.         go to 1
  29.     2   continue
  30. * on va reduire le maillage
  31.         ipt1=imamod
  32.         segact ipt1
  33.         segini,ipt2=ipt1
  34.         igard=0
  35.         nbelem=ipt1.num(/2)
  36.         if(idim.eq.3) then
  37. *
  38. *  boucle sur le maillage support en 3D
  39. *
  40.  
  41.           do 10 iel=1,nbelem
  42.             ip1=ipt1.num(2,iel)
  43.             ip2=ipt1.num(3,iel)
  44.             ip3=ipt1.num(4,iel)
  45.             xp1=xcoor((ip1-1)*(idim+1)+1)
  46.             yp1=xcoor((ip1-1)*(idim+1)+2)
  47.             zp1=xcoor((ip1-1)*(idim+1)+3)
  48.             xp2=xcoor((ip2-1)*(idim+1)+1)
  49.             yp2=xcoor((ip2-1)*(idim+1)+2)
  50.             zp2=xcoor((ip2-1)*(idim+1)+3)
  51.             xp3=xcoor((ip3-1)*(idim+1)+1)
  52.             yp3=xcoor((ip3-1)*(idim+1)+2)
  53.             zp3=xcoor((ip3-1)*(idim+1)+3)
  54. *
  55. *  verification si autre point dans la zone de selection
  56.             jp=ipt1.num(5,iel)
  57. *      write (6,*) 'mail ',ip1,ip2,ip3,jp
  58.             xp =xcoor((jp-1)*(idim+1)+1)
  59.             yp =xcoor((jp-1)*(idim+1)+2)
  60.             zp =xcoor((jp-1)*(idim+1)+3)
  61. *  verification que pas relation du point sur lui meme
  62.             if (jp.eq.ip1) goto 10
  63.             if (jp.eq.ip2) goto 10
  64.             if (jp.eq.ip3) goto 10
  65. *  verif distance au centre de gravite
  66.             xg=(xp1+xp2+xp3)/3d0
  67.             yg=(yp1+yp2+yp3)/3d0
  68.             zg=(zp1+zp2+zp3)/3d0
  69.             d1=sqrt((xg-xp1)**2+(yg-yp1)**2+(zg-zp1)**2)
  70.             d2=sqrt((xg-xp2)**2+(yg-yp2)**2+(zg-zp2)**2)
  71.             d3=sqrt((xg-xp3)**2+(yg-yp3)**2+(zg-zp3)**2)
  72.             dp=sqrt((xg-xp)**2+(yg-yp)**2+(zg-zp)**2)/2d0
  73.             if (dp.gt.d1.and.dp.gt.d2.and.dp.gt.d3) goto 10
  74. *       write (6,*) 'contact test distance ok',dp,d1,d2,d3
  75. *  verif position par rapport aux cotés
  76. *  cote 1 2
  77.             xv=xp2-xp1
  78.             yv=yp2-yp1
  79.             zv=zp2-zp1
  80.             dv=sqrt((xv**2)+(yv**2)+(zv**2))
  81.             if(.not.(dv.lt.1d-10).and.
  82.      >         .not.(dv.gt.1d-10))
  83.      >        write (6,*) 'prob 1 dans impo32'
  84.             if (abs(dv).le.xpetit) write (6,*) 'prob 1 imp32'
  85.             xv=xv/dv
  86.             yv=yv/dv
  87.             zv=zv/dv
  88.             scal=(xp-xp1)*xv+(yp-yp1)*yv+(zp-zp1)*zv
  89.             xm=xp1+scal*xv
  90.             ym=yp1+scal*yv
  91.             zm=zp1+scal*zv
  92.             scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
  93.             xn=(xg-xm)-scal*xv
  94.             yn=(yg-ym)-scal*yv
  95.             zn=(zg-zm)-scal*zv
  96.             dn=sqrt(xn**2+yn**2+zn**2)
  97.             scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  98.             dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  99.             if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
  100. *       write (6,*) 'contact test position 1 ok',scal/(dn*dpm),scal,
  101. *    >  dn,dpm
  102. *  cote 2 3
  103.             xv=xp3-xp2
  104.             yv=yp3-yp2
  105.             zv=zp3-zp2
  106.             dv=sqrt((xv**2)+(yv**2)+(zv**2))
  107.             if(.not.(dv.lt.1d-10).and.
  108.      >         .not.(dv.gt.1d-10))
  109.      >         write (6,*) 'prob 2 dans impo32'
  110.             if (abs(dv).le.xpetit) write (6,*) 'prob 2 imp32'
  111.             xv=xv/dv
  112.             yv=yv/dv
  113.             zv=zv/dv
  114.             scal=(xp-xp2)*xv+(yp-yp2)*yv+(zp-zp2)*zv
  115.             xm=xp2+scal*xv
  116.             ym=yp2+scal*yv
  117.             zm=zp2+scal*zv
  118.             scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
  119.             xn=(xg-xm)-scal*xv
  120.             yn=(yg-ym)-scal*yv
  121.             zn=(zg-zm)-scal*zv
  122.             dn=sqrt(xn**2+yn**2+zn**2)
  123.             scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  124.             dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  125.             if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
  126. *          write (6,*) 'contact test position 2 ok',scal/(dn*dpm),scal,
  127. *    >  dn,dpm
  128. *  cote 3 1
  129.             xv=xp1-xp3
  130.             yv=yp1-yp3
  131.             zv=zp1-zp3
  132.             dv=sqrt((xv**2)+(yv**2)+(zv**2))
  133.             if(.not.(dv.lt.1d-10).and.
  134.      >         .not.(dv.gt.1d-10))
  135.      >         write (6,*) 'prob 3 dans impo32'
  136.             if (abs(dv).le.xpetit) write (6,*) 'prob 3 imp32'
  137.             xv=xv/dv
  138.             yv=yv/dv
  139.             zv=zv/dv
  140.             scal=(xp-xp3)*xv+(yp-yp3)*yv+(zp-zp3)*zv
  141.             xm=xp3+scal*xv
  142.             ym=yp3+scal*yv
  143.             zm=zp3+scal*zv
  144.             scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
  145.             xn=(xg-xm)-scal*xv
  146.             yn=(yg-ym)-scal*yv
  147.             zn=(zg-zm)-scal*zv
  148.             dn=sqrt(xn**2+yn**2+zn**2)
  149.             scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  150.             dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  151.             if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
  152. *       write (6,*) 'contact test position 3 ok',scal/(dn*dpm),scal,
  153. *    >  dn,dpm
  154. * on a un point ou imposer la relation
  155. * Quelle est la relation ???
  156. *     write (6,*) 'contact potentiel '
  157. *
  158. *   direction de la relation
  159. *
  160.  
  161.             xn=(yp1-yp2)*(zp2-zp3)-(zp1-zp2)*(yp2-yp3)
  162.             yn=(zp1-zp2)*(xp2-xp3)-(xp1-xp2)*(zp2-zp3)
  163.             zn=(xp1-xp2)*(yp2-yp3)-(yp1-yp2)*(xp2-xp3)
  164.             dn=sqrt(xn**2+yn**2+zn**2)
  165.             if(.not.(dn.lt.1d-10).and.
  166.      >         .not.(dn.gt.1d-10))
  167.      >         write (6,*) 'prob 4 dans impo32'
  168.             if (abs(dn).le.xpetit) write (6,*) 'prob 4 imp32'
  169.             xn=xn/dn
  170.             yn=yn/dn
  171.             zn=zn/dn
  172.             if (abs(xn).lt.1d-10) xn=0.d0
  173.             if (abs(yn).lt.1d-10) yn=0.d0
  174.             if (abs(zn).lt.1d-10) zn=0.d0
  175. *
  176. *   recherche du point a mettre en relation (coordonnees barycentriques)
  177. *
  178.             scal=(xp-xg)*xn+(yp-yg)*yn+(zp-zg)*zn
  179. *
  180. *  verif bon cote (a la tolerance pres)
  181. *
  182.             if (scal.lt.-0.2d0*d1.and.scal.lt.-0.2d0*d2.and.
  183.      >          scal.lt.-0.2d0*d3) goto 10
  184.  
  185. *
  186. *   on range dans le meleme
  187. *
  188.  
  189.             igard=igard+1
  190.             ipt2.num(1,igard)=ipt1.num(1,iel)
  191.             ipt2.num(2,igard)=ip1
  192.             ipt2.num(3,igard)=ip2
  193.             ipt2.num(4,igard)=ip3
  194.             ipt2.num(5,igard)=jp
  195.    10     continue
  196.           segdes imodel,ipt1
  197.           if(igard.eq.nbelem) then
  198.             mmode2.kmodel(imo-imoin)=imodel
  199.             segsup ipt2
  200.             go to 1
  201.           endif
  202.           if( igard.eq.0) then
  203.             imoin=imoin+1
  204.             segsup ipt2
  205.             go to 1
  206.           endif
  207.           nbelem=igard
  208.           nbnn=ipt2.num(/1)
  209.           nbref=0
  210.           nbsous=0
  211.           segadj ipt2
  212.           segini,imode2=imodel
  213.           imode2.imamod=ipt2
  214.           kmodel(imo-imoin)=imode2
  215.           segdes imode2,ipt2
  216.           go to 1
  217.         elseif(idim.eq.2) then
  218.           if( ipt1.num(/1).eq.5) then
  219. *  formulation faible
  220.             do 210 iel=1,ipt1.num(/2)
  221.               ip1=ipt1.num(2,iel)
  222.               ip2=ipt1.num(3,iel)
  223.               ip3=ipt1.num(4,iel)
  224.               ip4=ipt1.num(5,iel)
  225.               xp1=xcoor((ip1-1)*(idim+1)+1)
  226.               yp1=xcoor((ip1-1)*(idim+1)+2)
  227.               xp2=xcoor((ip2-1)*(idim+1)+1)
  228.               yp2=xcoor((ip2-1)*(idim+1)+2)
  229.               xp3=xcoor((ip3-1)*(idim+1)+1)
  230.               yp3=xcoor((ip3-1)*(idim+1)+2)
  231.               xp4=xcoor((ip4-1)*(idim+1)+1)
  232.               yp4=xcoor((ip4-1)*(idim+1)+2)
  233. *
  234. *   critere d'acceptation de l'éléments:
  235. *   une extrémité d'un segment dans le cercle de sélection de l'autre segment
  236. *   orientations respectives correctes des 2 segments
  237. *
  238. *  verification (provisoire) que pas de noeuds doubles.
  239.               if (ip1.eq.ip3) goto 210
  240.               if (ip1.eq.ip4) goto 210
  241.               if (ip2.eq.ip3) goto 210
  242.               if (ip2.eq.ip4) goto 210
  243. *
  244.               xl13=sqrt((xp3-xp1)**2+(yp3-yp1)**2)
  245.               xl14=sqrt((xp4-xp1)**2+(yp4-yp1)**2)
  246.               xl23=sqrt((xp3-xp2)**2+(yp3-yp2)**2)
  247.               xl24=sqrt((xp4-xp2)**2+(yp4-yp2)**2)
  248.               sca312=(xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2)
  249.               sca412=(xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2)
  250.               sca134=(xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4)
  251.               sca234=(xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4)
  252.               if (sca312/(xl13*xl23).gt.0.7.and.
  253.      >            sca412/(xl14*xl24).gt.0.7.and.
  254.      >            sca134/(xl13*xl14).gt.0.7.and.
  255.      >            sca234/(xl23*xl24).gt.0.7)   goto 210
  256. *
  257.               scal=(xp1-xp2)*(xp3-xp4)+(yp1-yp2)*(yp3-yp4)
  258.               if (scal.gt.0.) goto 210
  259.  
  260. * on a un element ou imposer la relation
  261. * Quelle est la relation ???
  262. *
  263. *   direction du contact
  264. *
  265.               xr=(xp2-xp1+xp3-xp4)
  266.               yr=(yp2-yp1+yp3-yp4)
  267.               sr2 = xr**2 + yr**2
  268.               sr = sqrt (sr2)
  269.               xr=xr/sr
  270.               yr=yr/sr
  271. *
  272. *   projection des points sur la direction du contact
  273. *
  274.               xl1=xp1*xr+yp1*yr
  275.               xl2=xp2*xr+yp2*yr
  276.               xl3=xp3*xr+yp3*yr
  277.               xl4=xp4*xr+yp4*yr
  278. *
  279. *   distance des points à leurs projections
  280. *
  281.               d1=(xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr
  282.               d2=(xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr
  283.               d3=(xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr
  284.               d4=(xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr
  285. *
  286. *   calcul de l'intersection des projections.
  287. *
  288.               xm1=min(xl1,xl2)
  289.               xm2=max(xl1,xl2)
  290.               xm3=min(xl3,xl4)
  291.               xm4=max(xl3,xl4)
  292. *
  293. *  critere
  294. *
  295.               xcr=max(xm2-xm1,xm4-xm3)*1d-4
  296.               xi=max(xm1,xm3)
  297.               xj=min(xm2,xm4)
  298. *
  299.               if (xi.ge.xj-xcr) goto 210
  300.  
  301. *   on range dans le meleme
  302. *
  303.               igard=igard+1
  304.               ipt2.num(1,igard)=ipt1.num(1,iel)
  305.               ipt2.num(2,igard)=ip1
  306.               ipt2.num(3,igard)=ip2
  307.               ipt2.num(4,igard)=ip3
  308.               ipt2.num(5,igard)=ip4
  309. *
  310.   210       continue
  311.             segdes imodel,ipt1
  312.             if(igard.eq.nbelem) then
  313.               mmode2.kmodel(imo-imoin)=imodel
  314.               segsup ipt2
  315.               go to 1
  316.             endif
  317.             if( igard.eq.0) then
  318.               imoin=imoin+1
  319.               segsup ipt2
  320.               go to 1
  321.             endif
  322.             nbelem=igard
  323.             nbnn=ipt2.num(/1)
  324.             nbref=0
  325.             nbsous=0
  326.             segadj ipt2
  327.             segini,imode2=imodel
  328.             imode2.imamod=ipt2
  329.             kmodel(imo-imoin)=imode2
  330.             segdes imode2,ipt2
  331.             go to 1
  332.           elseif(ipt1.num(/1).eq.4) then
  333. * formulation 1 point - 1 segment
  334.             do 110 iel=1,ipt1.num(/2)
  335.               ip1=ipt1.num(2,iel)
  336.               ip2=ipt1.num(3,iel)
  337.               xp1=xcoor((ip1-1)*(idim+1)+1)
  338.               yp1=xcoor((ip1-1)*(idim+1)+2)
  339.               xp2=xcoor((ip2-1)*(idim+1)+1)
  340.               yp2=xcoor((ip2-1)*(idim+1)+2)
  341. *
  342. *  verification si autre point dans le cercle de selection
  343.               jp=ipt1.num(4,iel)
  344. *  verification que pas relation du point sur lui meme
  345.               if (jp.eq.ip1) goto 110
  346.               if (jp.eq.ip2) goto 110
  347.               xp=xcoor((jp-1)*(idim+1)+1)
  348.               yp=xcoor((jp-1)*(idim+1)+2)
  349.               scal=(xp-xp1)*(xp-xp2)+(yp-yp1)*(yp-yp2)
  350.               xl1=sqrt((xp-xp1)**2+(yp-yp1)**2)
  351.               xl2=sqrt((xp-xp2)**2+(yp-yp2)**2)
  352.               scal=scal/(xl1*xl2)
  353.               if (scal.gt.0.7) goto 110
  354.               igard=igard+1
  355.               ipt2.num(1,igard)=ipt1.num(1,iel)
  356.               ipt2.num(2,igard)=ip1
  357.               ipt2.num(3,igard)=ip2
  358.               ipt2.num(4,igard)=jp
  359.   110       continue
  360.             segdes imodel,ipt1
  361.             if(igard.eq.nbelem) then
  362.               mmode2.kmodel(imo-imoin)=imodel
  363.               segsup ipt2
  364.               go to 1
  365.             endif
  366.             if( igard.eq.0) then
  367.               imoin=imoin+1
  368.               segsup ipt2
  369.               go to 1
  370.             endif
  371.             nbelem=igard
  372.             nbnn=ipt2.num(/1)
  373.             nbref=0
  374.             nbsous=0
  375.             segadj ipt2
  376.             segini,imode2=imodel
  377.             imode2.imamod=ipt2
  378.             mmode2.kmodel(imo-imoin)=imode2
  379.             segdes imode2,ipt2
  380.             go to 1
  381.           endif
  382.         endif
  383.     1 continue
  384.       segdes mmodel
  385.       if(imoin.ne.0) then
  386.          n1=kmodel(/1)-imoin
  387.          if(n1.eq.0) then
  388.             call erreur (21)
  389.             return
  390.          endif
  391.          segadj mmode2
  392.       endif
  393.       segdes mmode2
  394.       return
  395.       end
  396.  
  397.  
  398.  

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