Télécharger impos2.eso

Retour à la liste

Numérotation des lignes :

  1. C IMPOS2 SOURCE PV 16/07/20 21:15:02 9040
  2.  
  3. * impo bloc en 2D
  4.  
  5. SUBROUTINE IMPOS2
  6.  
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10. -INC CCOPTIO
  11.  
  12. -INC SMCOORD
  13. -INC SMELEME
  14. -INC SMRIGID
  15. -INC SMCHPOI
  16. -INC SMCHAML
  17. -INC CCREEL
  18.  
  19. * directions moyennes aux sommets
  20. segment mfopa2
  21. real*8 xms(nbpts),yms(nbpts)
  22. endsegment
  23.  
  24. character*4 modepl(4),moforc(4)
  25.  
  26. data modepl /'UX ','UY ','UR ','UZ '/
  27. data moforc /'FX ','FY ','FR ','FZ '/
  28. *
  29. idimp1 = IDIM + 1
  30. *
  31. * Lecture du maillage support des conditions de contact
  32. *
  33. call lirobj('MAILLAGE',ipt1,1,iretou)
  34. call lirent(isous,0,iretou)
  35. mchel1=0
  36. call lirent(mchel1,0,iretou)
  37. if (ierr.ne.0) return
  38. *
  39. * Activation du maillage et petites verifications
  40. *
  41. segact ipt1
  42. nbno1 = ipt1.num(/1)
  43. nbel1 = ipt1.num(/2)
  44. if (ipt1.lisous(/1).ne.0) call erreur(25)
  45. if (ipt1.itypel.ne.22) call erreur(16)
  46. if (nbno1.ne.4 .and. nbno1.ne.5) call erreur(16)
  47. if (ierr.ne.0) goto 900
  48. *
  49. * Indice des ddls concernes en fonction du mode 2D
  50. *
  51. imo = 1
  52. if (ifour.eq.0) imo = 3
  53. *
  54. * Increment sur le nombre d'elements de contact retenus et utilises
  55. * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les
  56. * longueurs des segments adequats
  57. incnel = 2000
  58. *
  59. * Creation de la raideur des conditions de contact
  60. * Remplissage de l'entete commun
  61. *
  62. nrigel = 1
  63. segini mrigid
  64. ichole = 0
  65. imgeo1 = 0
  66. imgeo2 = 0
  67. isupeq = 0
  68. iforig = ifomod
  69. coerig(1)=1.d0
  70. mtymat='RIGIDITE'
  71. *
  72. * Creation du chpoint de depi
  73. *
  74. nat=1
  75. nsoupo=1
  76. segini mchpoi
  77. mtypoi='DEPIMP'
  78. mochde='engendré par impose'
  79. ifopoi=ifomod
  80. jattri(1)=2
  81. *
  82. nc=1
  83. segini msoupo
  84. ipchp(1)=msoupo
  85. nocomp(1)='FLX '
  86. nocons(1)=' '
  87. *
  88. nbnn =1
  89. nbelem=incnel
  90. nbref =0
  91. nbsous=0
  92. segini ipt7
  93. ipt7.itypel=1
  94. igeoc=ipt7
  95. *
  96. n=incnel
  97. segini mpoval
  98. ipoval=mpoval
  99. * Calcul des directions moyennes
  100. segact mcoord
  101. nbpts = xcoor(/1)/(idim+1)
  102. segini mfopa2
  103. ip1o=0
  104. ip2o=0
  105. * write (6,*) ' impos2 iel1 ',nbel1,idimp1
  106. do 820 iel=1,nbel1
  107. *
  108. ip1=ipt1.num(2,iel)
  109. ip2=ipt1.num(3,iel)
  110. * on suppose que les elements de frottements appuyes sur la même facette se suivent
  111. * et on ne les comptabilise qu'une fois
  112. if (ip1.eq.ip1o.and.ip2.eq.ip2o) goto 820
  113. ip1o=ip1
  114. ip2o=ip2
  115. ipv1 = (ip1-1)*idimp1
  116. xp1 = xcoor(ipv1+1)
  117. yp1 = xcoor(ipv1+2)
  118. ipv2 = (ip2-1)*idimp1
  119. xp2 = xcoor(ipv2+1)
  120. yp2 = xcoor(ipv2+2)
  121. *
  122. * normale a la droite (12)
  123. *
  124. x12 = xp2 - xp1
  125. y12 = yp2 - yp1
  126. xn = -y12
  127. yn = x12
  128. sn = sqrt (xn*xn + yn*yn)
  129. sn = max(xpetit,sn)
  130. xn = xn/sn
  131. yn = yn/sn
  132. xms(ip1)=xms(ip1)+xn
  133. yms(ip1)=yms(ip1)+yn
  134. xms(ip2)=xms(ip2)+xn
  135. yms(ip2)=yms(ip2)+yn
  136. 820 continue
  137. do 821 ip = 1,nbpts
  138. sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)
  139. sn = max(sqrt(abs(sn)),xpetit*1d10)
  140. xms(ip)=xms(ip)/sn
  141. yms(ip)=yms(ip)/sn
  142. 821 continue
  143. *
  144.  
  145. *
  146. * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval
  147. nelch = 0
  148. *
  149. *=======================================================================
  150. * Formulation "faible" du contact : relation segment-segment (type 0)
  151. *=======================================================================
  152. * Nouvelle formulation, une seule relation par element maitre
  153. * Il faut donc reunir les relations portant sur le meme multiplicateur
  154. ************************************************************************
  155. if (nbno1.ne.5) goto 500
  156. *
  157. * Creation du meleme associe a la relation
  158. *
  159. if (mchel1.ne.0) then
  160. segact mchel1
  161. mcham1 = mchel1.ichaml(isous)
  162. segact mcham1
  163. melva1=mcham1.ielval(1)
  164. segact melva1
  165. nptel = melva1.velche(/1)
  166. nel = melva1.velche(/2)
  167. nptel1=min(nptel,4)
  168. endif
  169.  
  170. nbnn =5
  171. nbelem=incnel
  172. nbsous=0
  173. nbref =0
  174. segini meleme
  175. itypel=22
  176. irigel(1,nrigel) = meleme
  177. *
  178. * Creation du descriptif commun a toutes les raideurs
  179. *
  180. nligrp=9
  181. nligrd=nligrp
  182. segini,descr
  183. lisinc(1)='LX '
  184. lisdua(1)='FLX '
  185. noelep(1)=1
  186. noeled(1)=1
  187. do 100 i = 2, nligrp, 2
  188. lisinc(i )=modepl(imo)
  189. lisinc(i+1)=modepl(imo+1)
  190. lisdua(i )=moforc(imo)
  191. lisdua(i+1)=moforc(imo+1)
  192. noelep(i )=(i+2)/2
  193. noelep(i+1)=noelep(i)
  194. noeled(i )=noelep(i)
  195. noeled(i+1)=noelep(i)
  196. 100 continue
  197. segdes,descr
  198. irigel(3,nrigel) = descr
  199. *
  200. * creation du segment xmatri
  201. *
  202. nelrig=incnel
  203. segini xmatri
  204. irigel(4,nrigel) = xmatri
  205. *
  206. * ce qu'on cree est unilateral
  207. *
  208. irigel(6,nrigel) = 1
  209. *
  210. * ce qu'on cree est symetrique
  211. *
  212. irigel(7,nrigel) = 0
  213. *
  214. * Nombre d'elements crees dans meleme=irigel(nrigel,1), ipt7 et mpoval
  215. nelri0 = 0
  216. *
  217. * boucle sur le maillage support
  218. *
  219. do 110 iel = 1, nbel1
  220. *
  221. xjr = 0d0
  222. if (mchel1.ne.0) then
  223. nel1 = min (iel,nel)
  224. xjr = melva1.velche(nptel1,nel1)
  225. endif
  226. ip1 = ipt1.num(2,iel)
  227. ip2 = ipt1.num(3,iel)
  228. ip3 = ipt1.num(4,iel)
  229. ip4 = ipt1.num(5,iel)
  230. * write(ioimp,*) iel,ip1,ip2,ip3,ip4,ipt1.num(1,iel)
  231. * verification (provisoire) que pas de noeuds doubles
  232. if (ip1.eq.ip3) goto 110
  233. if (ip1.eq.ip4) goto 110
  234. if (ip2.eq.ip3) goto 110
  235. if (ip2.eq.ip4) goto 110
  236. *
  237. ipv = (ip1-1)*idimp1
  238. xp1 = xcoor(ipv+1)
  239. yp1 = xcoor(ipv+2)
  240. ipv = (ip2-1)*idimp1
  241. xp2 = xcoor(ipv+1)
  242. yp2 = xcoor(ipv+2)
  243. xp12 = xp2 - xp1
  244. yp12 = yp2 - yp1
  245. *
  246. ipv = (ip3-1)*idimp1
  247. xp3 = xcoor(ipv+1)
  248. yp3 = xcoor(ipv+2)
  249. ipv = (ip4-1)*idimp1
  250. xp4 = xcoor(ipv+1)
  251. yp4 = xcoor(ipv+2)
  252. xp34 = xp4 - xp3
  253. yp34 = yp4 - yp3
  254. *
  255. * orientations respectives correctes des 2 segments :
  256. * "normales de sens opposes" = produit scalaire negatif ou nul
  257. scal = xp12*xp34 + yp12*yp34
  258. xl12 = sqrt(xp12**2+yp12**2)
  259. xl34 = sqrt(xp34**2+yp34**2)
  260. * if (scal.gt.-xl12*xl34*0.5) goto 110
  261. if (scal.gt.0.d0) goto 110
  262. *
  263. * critere d'acceptation de l'élément :
  264. * angles des diagonales
  265. xl13 = sqrt((xp3-xp1)**2+(yp3-yp1)**2)
  266. xl14 = sqrt((xp4-xp1)**2+(yp4-yp1)**2)
  267. xl23 = sqrt((xp3-xp2)**2+(yp3-yp2)**2)
  268. xl24 = sqrt((xp4-xp2)**2+(yp4-yp2)**2)
  269. sca312 = (xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2)
  270. sca412 = (xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2)
  271. sca134 = (xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4)
  272. sca234 = (xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4)
  273. * write(ioimp,*) sca312/(xl13*xl23),sca412/(xl14*xl24),
  274. * & sca134/(xl13*xl14),sca234/(xl23*xl24)
  275. ** if (sca312/(xl13*xl23).gt.0.50.and.
  276. ** > sca412/(xl14*xl24).gt.0.50.and.
  277. ** > sca134/(xl13*xl14).gt.0.50.and.
  278. ** > sca234/(xl23*xl24).gt.0.50) goto 110
  279.  
  280. * nouveau critere acceptation
  281. * pts dans un cercle centree sur 1-2 aggrandi
  282. xp1e=xp1-(xp2-xp1)
  283. yp1e=yp1-(yp2-yp1)
  284. xp2e=xp2+(xp2-xp1)
  285. yp2e=yp2+(yp2-yp1)
  286. sca312 = (xp3-xp1e)*(xp3-xp2e)+(yp3-yp1e)*(yp3-yp2e)
  287. sca412 = (xp4-xp1e)*(xp4-xp2e)+(yp4-yp1e)*(yp4-yp2e)
  288. if (sca312.gt.0.d0.and.
  289. > sca412.gt.0.d0) goto 110
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299. *
  300. * Quelle est la relation ???
  301. *
  302. * direction du contact unitaire
  303. *
  304. xr = xp12 - xp34
  305. yr = yp12 - yp34
  306. sr = sqrt(xr*xr + yr*yr)
  307. xr = xr/sr
  308. yr = yr/sr
  309. * write(ioimp,*) 'direction contact',xr,yr,yr,-xr
  310. *
  311. * projection des points sur la direction du contact
  312. *
  313. xl1 = xp1*xr + yp1*yr
  314. xl2 = xp2*xr + yp2*yr
  315. xl3 = xp3*xr + yp3*yr
  316. xl4 = xp4*xr + yp4*yr
  317.  
  318. * write(ioimp,*) 'projection pts sur contact',xl1,xl2,xl3,xl4
  319. *
  320. * calcul de l'intersection des projections
  321. *
  322. xm1 = min(xl1,xl2)
  323. xm2 = max(xl1,xl2)
  324. xm3 = min(xl3,xl4)
  325. xm4 = max(xl3,xl4)
  326. * write(ioimp,*) ' xmi',xm1,xm2,xm3,xm4
  327. *
  328. * critere de precision sur l'intersection
  329. *
  330. xcr = min(xm2-xm1,xm4-xm3)*(1.d-10)
  331. *
  332. * taille de l'intersection
  333. *
  334. xi = max(xm1,xm3)
  335. xj = min(xm2,xm4)
  336. xl = xj - xi
  337. * write(ioimp,*) ' intersection',xi,xj,xl,xcr
  338. *
  339. * write (6,*) ' impos2 ',ip1,ip2,ip3,ip4,xcr,xl
  340. * if (xl.le.xcr) goto 110
  341. if (xl.le.0.d0) goto 110
  342. *
  343. * distance des points a leur projection
  344. *
  345. d1 = (xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr
  346. d2 = (xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr
  347. d3 = (xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr
  348. d4 = (xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr
  349. *
  350. * coordonnées paramétriques de l'intersection sur les segments 1-2 et 3-4
  351. *
  352. xi2 = (xi-xl1) / (xl2-xl1)
  353. xi1 = (xl2-xi) / (xl2-xl1)
  354.  
  355. xj2 = (xj-xl1) / (xl2-xl1)
  356. xj1 = (xl2-xj) / (xl2-xl1)
  357.  
  358. xi4 = (xi-xl3) / (xl3-xl4)
  359. xi3 = (xl4-xi) / (xl3-xl4)
  360.  
  361. xj4 = (xj-xl3) / (xl3-xl4)
  362. xj3 = (xl4-xj) / (xl3-xl4)
  363. *
  364. * write(ioimp,*) ' xi1 xi2 xi3 xi4 xj1 xj2 xj3 xj4 xl '
  365. * write(ioimp,*) xi1,xi2,xi3,xi4,xj1,xj2,xj3,xj4,xl
  366. * write(ioimp,*) ' d1 d2 d3 d4 ',d1,d2,d3,d4
  367. *
  368. * surface actuelle
  369. *
  370. sc = ((xi1+xj1)*d1+(xi2+xj2)*d2+(xi3+xj3)*d3+(xi4+xj4)*d4)*0.5
  371. * write(ioimp,*) 'Surface actuelle :',sc
  372. *
  373. * on a un element ou imposer la relation a ajouter
  374. *
  375. nelri0 = nelri0 + 1
  376. nelch = nelch +1
  377. *
  378. * on ajuste les differents segments
  379. *
  380. if (nelri0.gt.nelrig) then
  381. nelrig = nelrig + incnel
  382. segadj,xmatri
  383. nbelem = nbelem + incnel
  384. nbnn = 5
  385. segadj,meleme
  386. nbnn = 1
  387. segadj,ipt7
  388. n = n + incnel
  389. segadj,mpoval
  390. endif
  391. *
  392. * on range dans le meleme
  393. *
  394. num(1,nelri0) = ipt1.num(1,iel)
  395. num(2,nelri0) = ip1
  396. num(3,nelri0) = ip2
  397. num(4,nelri0) = ip3
  398. num(5,nelri0) = ip4
  399. icolor(nelri0)=1
  400. *
  401. * on remplit le xmatri
  402. *
  403. * lambda
  404. re(1,1,nelri0)= 0.d0
  405. * ip1
  406. re(2,1,nelri0)= yr * (xi1+xj1) * 0.5 * xl
  407. re(3,1,nelri0)= -xr * (xi1+xj1) * 0.5 * xl
  408. * ip2
  409. re(4,1,nelri0)= yr * (xi2+xj2) * 0.5 * xl
  410. re(5,1,nelri0)= -xr * (xi2+xj2) * 0.5 * xl
  411. * ip3
  412. re(6,1,nelri0)= yr * (xi3+xj3) * 0.5 * xl
  413. re(7,1,nelri0)= -xr * (xi3+xj3) * 0.5 * xl
  414. * ip4
  415. re(8,1,nelri0)= yr * (xi4+xj4) * 0.5 * xl
  416. re(9,1,nelri0)= -xr * (xi4+xj4) * 0.5 * xl
  417. *
  418. * on transpose
  419. do 120 ic = 2, nligrp
  420. re(1,ic,nelri0)=re(ic,1,nelri0)
  421. 120 continue
  422. *
  423. * Le reste est nul
  424. *
  425. * remplissage du champoint de depi et du maillage support
  426. *
  427. ipt7.num(1,nelch)=ipt1.num(1,iel)
  428. vpocha(nelch,1) = (-sc - xjr) * xl
  429. * write (6,*) 'impos2 vpocha re ',vpocha(nelch,1),
  430. * > (re(ip,1,nelri0),ip=2,9)
  431. *
  432. 110 continue
  433.  
  434. * write (ioimp,*) ' nb relation type 0 ',nbelem,n,nelri0
  435.  
  436. * Ajustement au plus juste puis desactivation des segments lies
  437. * a la rigidite du type 0
  438. if (nelri0.ne.nelrig) then
  439. nelrig = nelri0
  440. segadj,xmatri
  441. nbelem = nelri0
  442. nbnn = 5
  443. segadj,meleme
  444. endif
  445. segdes,meleme,xmatri
  446. *
  447. * S'il n'y a pas d'elements en contact, alors pas de relation unilaterale
  448. * if (nelri0.eq.0) irigel(6,nrigel)=0
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455. GOTO 1000
  456.  
  457. *=======================================================================
  458. * Formulation "forte" (standard) du contact :
  459. *=======================================================================
  460. 500 continue
  461. *
  462. if (mchel1.ne.0) then
  463. segact mchel1
  464. mcham1 = mchel1.ichaml(isous)
  465. segact mcham1
  466. melva1=mcham1.ielval(1)
  467. segact melva1
  468. nptel = melva1.velche(/1)
  469. nel = melva1.velche(/2)
  470. nptel1=min(nptel,4)
  471. endif
  472. * Relation type 2 : noeud-segment
  473. *---------------------------------
  474. * creation du meleme associe a la relation
  475. *
  476. nbnn = 4
  477. nbelem = incnel
  478. nbsous = 0
  479. nbref = 0
  480. segini meleme
  481. itypel=22
  482. irigel(1,nrigel)=meleme
  483. *
  484. * creation du descriptif commun a toutes les raideurs
  485. *
  486. nligrp=7
  487. nligrd=nligrp
  488. segini descr
  489. lisinc(1)='LX '
  490. lisdua(1)='FLX '
  491. noelep(1)=1
  492. noeled(1)=1
  493. do 510 i=2,nligrp,2
  494. lisinc(i )=modepl(imo)
  495. lisinc(i+1)=modepl(imo+1)
  496. lisdua(i )=moforc(imo)
  497. lisdua(i+1)=moforc(imo+1)
  498. noelep(i )=(i+2)/2
  499. noelep(i+1)=noelep(i)
  500. noeled(i )=noelep(i)
  501. noeled(i+1)=noelep(i)
  502. 510 continue
  503. segdes,descr
  504. irigel(3,nrigel)=descr
  505. *
  506. * creation du segment xmatri
  507. *
  508. nelrig=incnel
  509. segini xmatri
  510. irigel(4,nrigel)=xmatri
  511. *
  512. * ce qu'on cree est unilateral
  513. *
  514. irigel(6,nrigel)=1
  515. *
  516. * ce qu'on cree est symetrique
  517. *
  518. irigel(7,nrigel)=0
  519. *
  520. * Nombre d'elements dans la rigidite de type 2
  521. nelri2=0
  522. *
  523. * boucle sur le maillage support
  524. *
  525. do 520 iel = 1,nbel1
  526. xjr = 0.d0
  527. if (mchel1.ne.0) then
  528. nel1 = min (iel,nel)
  529. xjr = melva1.velche(nptel1,nel1)
  530. * write(ioimp,*) iel,xjr,mchel1
  531. endif
  532. ip1 = ipt1.num(2,iel)
  533. ip2 = ipt1.num(3,iel)
  534. jp = ipt1.num(4,iel)
  535. *D write(ioimp,*) iel,ip1,ip2,jp
  536. *
  537. * verification que pas relation du point sur lui meme
  538. if (jp.eq.ip1) goto 520
  539. if (jp.eq.ip2) goto 520
  540.  
  541. ipv = (ip1-1)*idimp1
  542. xp1 = xcoor(ipv+1)
  543. yp1 = xcoor(ipv+2)
  544. ipv = (ip2-1)*idimp1
  545. xp2 = xcoor(ipv+1)
  546. yp2 = xcoor(ipv+2)
  547. ipv = (jp-1)*idimp1
  548. xp = xcoor(ipv+1)
  549. yp = xcoor(ipv+2)
  550.  
  551. x12 = xp2 - xp1
  552. y12 = yp2 - yp1
  553. x1p = xp - xp1
  554. y1p = yp - yp1
  555. x2p = xp - xp2
  556. y2p = yp - yp2
  557.  
  558. * distance signee de p a la ligne 1-2
  559. dp12s = x1p * y12 - y1p * x12
  560. * write(ioimp,*) 'dist. signee',dp12s
  561. * verif que le point est du bon cote du segment
  562. * if (dp12s.lt.0.) goto 520
  563.  
  564. * verification si autre point dans le cercle de selection
  565. d1p = x1p*x1p + y1p*y1p
  566. d2p = x2p*x2p + y2p*y2p
  567. if (abs(d1p).lt.XPETIT) d1p=1
  568. if (abs(d2p).lt.XPETIT) d2p=1
  569. scal = (x1p*x2p + y1p*y2p) / sqrt(d1p*d2p)
  570. * write(ioimp,*) 'cos(angle_1p2)',scal
  571. if (scal.gt.0.7) goto 520
  572. *
  573. * on a un point ou imposer la relation
  574. * Quelle est la relation ???
  575. *
  576. * direction de la relation
  577. *
  578. * calcul du pt a mettre en relation avec ip : alpha ip2 + (1-alpha) ip1
  579. *
  580. sr2 = x12*x12 + y12*y12
  581. alpha = (x1p*x12 + y1p*y12) / sr2
  582. *
  583. sr = sqrt(sr2)
  584. xr = x12 / sr
  585. yr = y12 / sr
  586. * normale normalisee
  587. xn=(1.-alpha)*xms(ip1)+alpha*xms(ip2)
  588. yn=(1.-alpha)*yms(ip1)+alpha*yms(ip2)
  589. sn=sqrt(xn*xn+yn*yn)
  590. if (abs(sn).lt.xpetit*1d10) then
  591. * si normale nulle on utilise la direction du segment pour la recalculer
  592. xn = -y12
  593. yn = x12
  594. sn = sqrt (xn*xn + yn*yn)
  595. sn = max(xpetit,sn)
  596. endif
  597. xn = xn/sn
  598. yn = yn/sn
  599. 1954 continue
  600. *
  601. * on a un element ou imposer la relation a ajouter
  602. *
  603. nelri2 = nelri2 + 1
  604. nelch = nelch + 1
  605. *
  606. * on ajuste les differents segments
  607. *
  608. if (nelri2.gt.nelrig) then
  609. nelrig = nelrig + incnel
  610. segadj,xmatri
  611. nbelem = nbelem + incnel
  612. nbnn = 4
  613. segadj,meleme
  614. nbnn = 1
  615. segadj,ipt7
  616. n = n + incnel
  617. segadj,mpoval
  618. endif
  619. *
  620. * on range dans le meleme
  621. *
  622. num(1,nelri2)=ipt1.num(1,iel)
  623. num(2,nelri2)=ip1
  624. num(3,nelri2)=ip2
  625. num(4,nelri2)=jp
  626. *
  627. * on remplit le xmatri
  628. *
  629. * lambda
  630. re(1,1,nelri2)= 0.d0
  631. * ip1
  632. re(2,1,nelri2)= -xn * (1.-alpha) * sr
  633. re(3,1,nelri2)= -yn * (1.-alpha) * sr
  634. * ip2
  635. re(4,1,nelri2)= -xn * alpha * sr
  636. re(5,1,nelri2)= -yn * alpha * sr
  637. * jp
  638. re(6,1,nelri2)= xn * sr
  639. re(7,1,nelri2)= yn * sr
  640. * on transpose
  641. re(1,2,nelri2) = re(2,1,nelri2)
  642. re(1,3,nelri2) = re(3,1,nelri2)
  643. re(1,4,nelri2) = re(4,1,nelri2)
  644. re(1,5,nelri2) = re(5,1,nelri2)
  645. re(1,6,nelri2) = re(6,1,nelri2)
  646. re(1,7,nelri2) = re(7,1,nelri2)
  647. * le reste est nul
  648. *
  649. * remplissage du champoint de depi
  650. *
  651. xjeu1 = x1p*xn + y1p*yn
  652. xjeu2 = x2p*xn + y2p*yn
  653. xjeu = (1.-alpha)*xjeu1 + alpha * xjeu2
  654. ipt7.num(1,nelch) = ipt1.num(1,iel)
  655. vpocha(nelch,1) = (-xjeu - xjr)*sr
  656. * write(ioimp,*) ' jeu type 2 ',nelri2,nelch,vpocha(nelch,1),
  657. * & ipt7.num(1,nelch),ip1,ip2,jp
  658. 520 continue
  659. *
  660. *D write (ioimp,*) ' nb relation type 2 ',nelri2
  661. *
  662. * Ajustement au plus juste puis desactivation des segments lies
  663. * a la rigidite du type 2
  664. if (nelri2.ne.nelrig) then
  665. nelrig = nelri2
  666. segadj,xmatri
  667. nbelem = nelri2
  668. nbnn = 4
  669. segadj,meleme
  670. endif
  671. segdes,meleme,xmatri
  672. *
  673. * si il n'y a rien on dit que pas unilateral pour pas passer dans unilater
  674. * if (nbelem.eq.0) irigel(6,nrigel) = 0
  675. *
  676. 700 CONTINUE
  677. *
  678. * GOTO 1000
  679. *
  680. *----------------------------
  681. * on renvoie le resultat
  682. *----------------------------
  683. 1000 CONTINUE
  684. segsup mfopa2
  685. *
  686. * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7
  687. * puis desactivation du chpoint
  688. if (n.ne.nelch) then
  689. n = nelch
  690. segadj,mpoval
  691. nbnn = 1
  692. nbelem = nelch
  693. nbsous = 0
  694. nbref = 0
  695. segadj,ipt7
  696. endif
  697. segdes,mpoval,ipt7
  698. segdes,msoupo,mchpoi
  699. * Desctivation de la matrice de raideur de contact
  700. segdes,mrigid
  701. *
  702. * Reunion des relations portant sur le meme multiplicateur de lagrange
  703. *
  704. call impofu(MRIGID,MCHPOI)
  705. *
  706. * Ecriture des objets resultats
  707. call ecrobj('RIGIDITE',mrigid)
  708. call ecrobj('CHPOINT',mchpoi)
  709.  
  710. 900 continue
  711. segdes,ipt1
  712.  
  713. return
  714. end
  715.  
  716.  
  717.  
  718.  
  719.  
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  

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