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    CB215821  25/07/02    21:15:02     12306          
  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.          segini xmatri
  172.          irigel(4,1)=xmatri
  173. *
  174. *       ce qu'on cree est une égalité
  175. *
  176.          irigel(6,1)=0
  177. *
  178. *       ce qu'on cree est symetrique
  179. *
  180.          irigel(7,1)=0
  181. *
  182. *       boucle sur le maillage support
  183.          do 10 iel=1,ipt1.num(/2)
  184. *
  185. *        noeud maitre
  186. *
  187.             ip1=ipt1.num(2,iel)
  188.             xp1=xcoor((ip1-1)*(idim+1)+1)
  189.             yp1=xcoor((ip1-1)*(idim+1)+2)
  190. *
  191. *        deuxieme noeud maitre
  192. *
  193.             ip2=ipt1.num(3,iel)
  194.             xp2=xcoor((ip2-1)*(idim+1)+1)
  195.             yp2=xcoor((ip2-1)*(idim+1)+2)
  196.             xr=xp2-xp1
  197.             yr=yp2-yp1
  198. *
  199. *        noeud esclave
  200. *
  201.             jp=ipt1.num(4,iel)
  202.             xp=xcoor((jp-1)*(idim+1)+1)
  203.             yp=xcoor((jp-1)*(idim+1)+2)
  204. *
  205. *        test de recherche du contact
  206. *
  207.             if(xp2.gt.xp1) then
  208.                if(xp.lt.xp1 .or. xp.gt.xp2) goto 10
  209.             else
  210.                if(xp.lt.xp2 .or. xp.gt.xp1) goto 10
  211.             endif
  212.             if((xp-xp1)*yr-(yp-yp1)*xr.gt.0.d0) goto 10
  213. *
  214. *        determination du noeud defenseur fictif
  215. *
  216.             sr2=xr**2+yr**2
  217.             xa=xp2-xp
  218.             ya=yp2-yp
  219.             alpha=(xa*xr+ya*yr)/sr2
  220.             if(abs(alpha).lt.1D-3) alpha=0.D0
  221.             if(alpha.gt.0.999D0) alpha=1.D0
  222.             if(alpha.lt.-0.999D0) alpha=-1.D0
  223.             sr=sqrt(sr2)
  224.             xr=xr/sr
  225.             yr=yr/sr
  226. *
  227. *        creation du xmatri
  228. *
  229. *            segini xmatri
  230. *            imatri=irigel(4,1)
  231.             nelri=nelri+1
  232.                         if(nelri.gt.nelrig) then
  233.                          nelrig=nelrig+50
  234.              segadj xmatri
  235.                         endif
  236. *            imattt(nelrig)=xmatri
  237. *
  238. *        multiplicateur
  239. *
  240.             re(1,1,nelri)= 0.d0
  241. *
  242. *        1er noeud maitre ip1
  243. *
  244.             re(2,1,nelri)=  yr * alpha
  245.             re(3,1,nelri)= -xr * alpha
  246. *
  247. *        2eme noeud maitre ip2
  248. *
  249.             re(4,1,nelri)=  yr * (1-alpha)
  250.             re(5,1,nelri)= -xr * (1-alpha)
  251. *
  252. *        noeud esclave jp
  253.             re(6,1,nelri)= -yr
  254.             re(7,1,nelri)=  xr
  255. *
  256. *        test de conditionnement
  257.             do 90 il=1,7
  258.                if(abs(re(il,1,nelri)).lt.1D-6) re(il,1,nelri)=0.D0
  259.                if(re(il,1,nelri).gt.(1.D0-1D-6)) re(il,1,nelri)=1.D0
  260.                if(re(il,1,nelri).lt.(-1.D0+1D-6)) re(il,1,nelri)=-1.D0
  261.  90         continue
  262. *
  263. *        transposition du xmatri C
  264.             do 30 il=1,1
  265.                do 301 ic=1,7
  266.                   re(il,ic,nelri)=re(ic,il,nelri)
  267.  301           continue
  268.  30         continue
  269. *
  270. *        modification du meleme
  271. *
  272.             nbelem=num(/2)+1
  273.             nbnn=4
  274.             nbsous=0
  275.             nbref=0
  276.             segadj meleme
  277.             num(1,nbelem)=ipt1.num(1,iel)
  278.             num(2,nbelem)=ip1
  279.             num(3,nbelem)=ip2
  280.             num(4,nbelem)=jp
  281. *
  282. *            segdes xmatri
  283. *
  284.  10      continue
  285.         if(nelri.ne.nelrig) then
  286.                   nelrig=nelri
  287.                   segadj xmatri
  288.                 endif
  289. *
  290.       else
  291. *
  292. *       creation des descriptifs communs a toutes les raideurs
  293. *
  294.          imo=1
  295.          if(ifour.eq.0) imo=3
  296.          nligrp=5
  297.          nligrd=nligrp
  298.          segini descr
  299.          irigel(3,1)=descr
  300.          lisinc(1)='LX'
  301.          lisdua(1)='FLX'
  302.          noelep(1)=1
  303.          noeled(1)=1
  304.          do 70 i=2,5,2
  305.             lisinc(i)=modepl(imo)
  306.             lisinc(i+1)=modepl(imo+1)
  307.             lisdua(i)=moforc(imo)
  308.             lisdua(i+1)=moforc(imo+1)
  309.             noelep(i)=(i+2)/2
  310.             noelep(i+1)=noelep(i)
  311.             noeled(i)=noelep(i)
  312.             noeled(i+1)=noelep(i)
  313.  70      continue
  314. *
  315. *       creation du segment imatri
  316. *
  317.          nelrig=50
  318.                  nelri=0
  319.          segini xmatri
  320.          irigel(4,1)=xmatri
  321. *
  322. *       ce qu'on cree est une égalité
  323. *
  324.          irigel(6,1)=0
  325. *
  326. *       ce qu'on cree est symetrique
  327. *
  328.          irigel(7,1)=0
  329. *
  330. *       boucle sur le maillage support
  331. *
  332.          do 50 iel=1,ipt1.num(/2)
  333. *
  334. *        noeud maitre
  335. *
  336.             ip1=ipt1.num(2,iel)
  337.             xp1=xcoor((ip1-1)*(idim+1)+1)
  338.             yp1=xcoor((ip1-1)*(idim+1)+2)
  339. *
  340. *        deuxieme noeud maitre fictif
  341. *
  342.             if(iel.lt.ipt1.num(/2)) then
  343.                ip2=ipt1.num(2,iel+1)
  344.                jcm=1
  345.             else
  346.                ip2=ipt1.num(2,iel-1)
  347.                jcm=-1
  348.             endif
  349.             xp2=xcoor((ip2-1)*(idim+1)+1)
  350.             yp2=xcoor((ip2-1)*(idim+1)+2)
  351.             xr=(xp2-xp1)*jcm
  352.             yr=(yp2-yp1)*jcm
  353. *
  354. *        noeud esclave
  355. *
  356.             jp=ipt1.num(3,iel)
  357.             xp=xcoor((jp-1)*(idim+1)+1)
  358.             yp=xcoor((jp-1)*(idim+1)+2)
  359. *
  360. *        verification que pas relation du point sur lui meme
  361. *
  362.             if (jp.eq.ip1) goto 50
  363. *
  364. *        detection du contact : point materiel
  365. *
  366.             if (iret.eq.0) then
  367.                if(((xp-xp1)*yr-(yp-yp1)*xr).gt.0.d0) goto 50
  368.                vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
  369.                vn1=vec(1)/vnorm
  370.                vn2=vec(2)/vnorm
  371.             endif
  372. *
  373. *        detection du contact : projectile a nez plat
  374. *
  375.             if (iret.eq.1) then
  376.                tv=vec(1)/vec(2)
  377.                aa=(xp-xp1)*tv
  378.                if(-(yp+aa-yp1)*xr.gt.0.d0) goto 50
  379.                vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
  380.                if(abs(xp-xp1).gt.abs(vlarg*vec(2)/vnorm)) goto 50
  381.                vn1=vec(1)/vnorm
  382.                vn2=vec(2)/vnorm
  383.             endif
  384. *
  385. *        detection du contact : projectile a nez hemispherique
  386. *
  387.             if (iret.eq.2) then
  388.                vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
  389.                if(abs(xp-xp1).gt.vrayo*abs(vec(2))/vnorm) goto 50
  390.                ypo=(vrayo/sqrt(1+(vec(1)*vec(1)/vec(2)*vec(2))))+yp
  391.                xpo=(ypo-yp)*vec(1)/vec(2)+xp
  392.                yy=ypo-sqrt(vrayo*vrayo-(xpo-xp1)*(xpo-xp1))
  393.                if(-(yy-yp1)*xr.gt.0.d0) goto 50
  394.                vnn=sqrt((xpo-xp1)*(xpo-xp1)+(ypo-yp1)*(ypo-yp1))
  395.                vn1=(xp1-xpo)/vnn
  396.                vn2=(yp1-ypo)/vnn
  397.             endif
  398. *
  399. *        detection du contact : projectile a nez conique
  400. *
  401.             if (iret.eq.3) then
  402.                ial=nint(sign(1.D0,(xp-xp1)))
  403.                vangl2=vangl-ial*atan(vec(1)/vec(2))
  404.                aa=abs(xp-xp1)*tan(vangl2)
  405.                if(-(yp+aa-yp1)*xr.gt.0.d0) goto 50
  406.                vnorm=sqrt(vec(1)*vec(1)+vec(2)*vec(2))
  407.                if(abs(xp-xp1).gt.vlarg*abs(vec(2))/vnorm) goto 50
  408.                v1=cos(vangl)*vec(1)+ial*sin(vangl)*vec(2)
  409.                v2=-ial*sin(vangl)*vec(1)+cos(vangl)*vec(2)
  410.                vnn=sqrt(v1*v1+v2*v2)
  411.                vn1=v1/vnn
  412.                vn2=v2/vnn
  413.             endif
  414. *
  415. *        on cree le xmatri C
  416. *
  417. *            segini xmatri
  418. *            imatri=irigel(4,1)
  419.             nelri=nelri+1
  420.                         if(nelri.gt.nelrig) then
  421.                          nelrig=nelrig+50
  422.              segadj xmatri
  423.                         endif
  424. *            imattt(nelrig)=xmatri
  425. *
  426. *        on remplit le xmatri
  427. *
  428. *        multiplicateur
  429.             re(1,1,nelri)= 0.d0
  430. *
  431. *        noeud maitre ip1
  432.             re(2,1,nelri)= vn1
  433.             re(3,1,nelri)=-vn2
  434. *
  435. *        noeud esclave jp
  436.             re(4,1,nelri)=-vn1
  437.             re(5,1,nelri)= vn2
  438. *
  439. *        test de conditionnement
  440.             do 95 il=1,5
  441.                if(abs(re(il,1,nelri)).lt.1D-6) re(il,1,nelri)=0.D0
  442.                if(re(il,1,nelri).gt.(1.D0-1D-6)) re(il,1,nelri)=1.D0
  443.                if(re(il,1,nelri).lt.(-1.D0+1D-6)) re(il,1,nelri)=-1.D0
  444.  95         continue
  445. *
  446. *        transposition du xmatri C
  447.             do 40 il=1,1
  448.                do 401 ic=1,5
  449.                   re(il,ic,nelri)=re(ic,il,nelri)
  450.  401           continue
  451.  40         continue
  452. *
  453. *        modification du meleme
  454. *
  455.             nbelem=num(/2)+1
  456.             nbnn=3
  457.             nbsous=0
  458.             nbref=0
  459.             segadj meleme
  460.             num(1,nbelem)=ipt1.num(1,iel)
  461.             num(2,nbelem)=ip1
  462.             num(3,nbelem)=jp
  463. *
  464. *            segdes xmatri
  465. *
  466.  50      continue
  467.          if(nelrig.ne.nelri) then
  468.                    nelrig=nelri
  469.                    segadj xmatri
  470.                  endif
  471. *
  472.       endif
  473.  
  474.       SEGDES,MCOORD
  475. *
  476. *     desactivation des segments
  477. *
  478. *      segdes xmatri,ipt1,mrigid,meleme,descr
  479. *
  480. *     renvoi du resultat
  481. *
  482.       if(nbelem.eq.0) then
  483.          irigel(6,1)=0
  484.          call ecrobj('ANNULE',mrigid)
  485.       else
  486.          call actobj('RIGIDITE',mrigid,1)
  487.          call ecrobj('RIGIDITE',mrigid)
  488.       endif
  489. *
  490.       end
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  

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