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

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