Télécharger impos2.eso

Retour à la liste

Numérotation des lignes :

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

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