Télécharger impos7.eso

Retour à la liste

Numérotation des lignes :

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

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