Télécharger la plaquette Cast3M


fr en

Télécharger impos7.eso

Retour à la liste

Numérotation des lignes : 

impos7

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

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