Télécharger impo32.eso

Retour à la liste

Numérotation des lignes :

  1. C IMPO32 SOURCE PV 20/04/26 21:15:01 10589
  2.  
  3. * impo bloc en 3D
  4.  
  5. SUBROUTINE IMPO32(ipt6)
  6.  
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10.  
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC CCREEL
  14. -INC CCGEOME
  15.  
  16. -INC SMCOORD
  17. -INC SMELEME
  18. -INC SMRIGID
  19. -INC SMCHPOI
  20. -INC SMCHAML
  21.  
  22. * Contact en formulation forte :
  23. * Segments utiles dans les cas dits pathologiques
  24. segment mfopa1
  25. integer lpoin(nfopa1,3)
  26. real*8 xpoin(nfopa1,3)
  27. real*8 zje1(nfopa1)
  28. endsegment
  29. segment mfopa2
  30. integer lsegm(nfopa2,4)
  31. real*8 xsegm(nfopa2,3)
  32. real*8 zje2(nfopa2)
  33. endsegment
  34. * directions moyennes aux sommets
  35. segment mfopa3
  36. real*8 xms(nbpts),yms(nbpts),zms(nbpts)
  37. endsegment
  38.  
  39. * Contact en formulation faible :
  40. * Segment utile pour le calcul de l'intersection des 2 triangles
  41. * projetes dans le plan de contact intermediaire
  42. segment mfaible
  43. real*8 cPT1C(3,3),cPT2C(3,3), cPT1A(4,2),cPT2A(4,2)
  44. real*8 cPIn0(6,2), cPIn(6,2), test(6)
  45. real*8 vT1A(3,2), vT2A(3,2)
  46. real*8 SuT1A, SuT2A, SuPIn
  47. *DBG-F real*8 xGIn,yGIn, b1T1,b2T1,b3T1, b1T2,b2T2,b3T2
  48. endsegment
  49.  
  50. PARAMETER ( X1s3=0.333333333333333333333333333333333333333333D0 )
  51.  
  52. character*4 modepl(3),moforc(3)
  53.  
  54. data modepl /'UX ','UY ','UZ '/
  55. data moforc /'FX ','FY ','FZ '/
  56. *
  57. idimp1 = IDIM + 1
  58. *
  59. * Lecture du maillage support des conditions de contact
  60. *
  61. ipt1 =0
  62. call lirobj('MAILLAGE',ipt1,1,iretou)
  63. call actobj('MAILLAGE',ipt1,1)
  64. isous=0
  65. call lirent(isous,0,iretou)
  66. mchel1=0
  67. call lirent(mchel1,0,iretou)
  68. C
  69. if (ierr.ne.0) return
  70. *
  71. * Activation du maillage et petites verifications
  72. *
  73. segact ipt1
  74. nbno1 = ipt1.num(/1)
  75. nbel1 = ipt1.num(/2)
  76. if (ipt1.lisous(/1).ne.0) call erreur(25)
  77. if (ipt1.itypel.ne.22) call erreur(16)
  78. if (nbno1.ne.5 .and. nbno1.ne.7) call erreur(16)
  79. if (ierr.ne.0) goto 900
  80. *
  81. * Increment sur le nombre d'elements de contact retenus et utilises
  82. * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les
  83. * longueurs des segments adequats
  84. incnel = 2000
  85. *
  86. * Creation de la raideur des conditions de contact
  87. * Remplissage de l'entete commun
  88. *
  89. nrigel = 1
  90. segini,mrigid
  91. ichole = 0
  92. imgeo1 = 0
  93. imgeo2 = 0
  94. isupeq = 0
  95. iforig = ifomod
  96. coerig(1) = 1.
  97. mtymat='RIGIDITE'
  98. *
  99. * MCHAML materiau associe au MMODEL
  100. MELVA1 = 0
  101. MELVA2 = 0
  102. IF (MCHEL1.NE.0) THEN
  103. SEGACT, MCHEL1
  104. * recherche rapide d'un maillage correspondant dans le mchaml
  105. DO 210 n2 = 1,mchel1.imache(/1)
  106. * write(ioimp,*) ' n2 imache(n2) ipt1 ',n2,mchel1.imache(n2),ipt1
  107. if (mchel1.imache(n2).eq.ipt1) goto 220
  108. 210 continue
  109.  
  110. goto 230
  111. 220 continue
  112. MCHAM1 = MCHEL1.ICHAML(n2)
  113. SEGACT, MCHAM1
  114. NBCO = MCHAM1.NOMCHE(/2)
  115. DO JCMP = 1,NBCO
  116. IF (MCHAM1.NOMCHE(JCMP).EQ.'JEU') THEN
  117. MELVA1 = MCHAM1.IELVAL(JCMP)
  118. SEGACT, MELVA1
  119. NELJ = MELVA1.VELCHE(/2)
  120. NPTELJ = min(MELVA1.VELCHE(/1),4)
  121. ENDIF
  122. IF (MCHAM1.NOMCHE(JCMP).EQ.'ADHE') THEN
  123. MELVA2 = MCHAM1.IELVAL(JCMP)
  124. SEGACT, MELVA2
  125. NELA = MELVA2.VELCHE(/2)
  126. NPTELA = min(MELVA2.VELCHE(/1),4)
  127. ENDIF
  128. ENDDO
  129. ENDIF
  130. 230 continue
  131. *
  132. * Creation du chpoint de depi
  133. *
  134. nat=1
  135. nsoupo=1
  136. segini mchpoi
  137. mtypoi='DEPIMP'
  138. mochde='engendré par impose'
  139. ifopoi=ifomod
  140. jattri(1)=2
  141. *
  142. nc=1
  143. IF (MELVA2.NE.0) nc=2
  144. segini msoupo
  145. ipchp(1)=msoupo
  146. nocomp(1)='FLX '
  147. nocons(1)=' '
  148. IF (MELVA2.NE.0) THEN
  149. nocomp(2)='FADH'
  150. nocons(2)=' '
  151. ENDIF
  152. *
  153. nbnn =1
  154. nbelem=incnel
  155. nbref =0
  156. nbsous=0
  157. segini ipt7
  158. ipt7.itypel=1
  159. igeoc=ipt7
  160. *
  161. n=incnel
  162. segini mpoval
  163. ipoval=mpoval
  164. *
  165. * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval
  166. nelch = 0
  167.  
  168. *=======================================================================
  169. * Formulation "forte" (standard) du contact :
  170. *=======================================================================
  171. * Element de contact 3D a 5 noeuds
  172. if (nbno1.eq.7) goto 500
  173.  
  174. * Calcul des directions moyennes
  175. segact mcoord
  176. segini mfopa3
  177. ip1o=0
  178. ip2o=0
  179. ip3o=0
  180. do 820 iel=1,nbel1
  181. *
  182. ip1=ipt1.num(2,iel)
  183. ip2=ipt1.num(3,iel)
  184. ip3=ipt1.num(4,iel)
  185. * on suppose que les elements de frottements appuyes sur la même facette se suivent
  186. * et on ne les comptabilise qu'une fois
  187. if (ip1.eq.ip1o.and.ip2.eq.ip2o.and.ip3.eq.ip3o) goto 820
  188. ip1o=ip1
  189. ip2o=ip2
  190. ip3o=ip3
  191. ipv1 = (ip1-1)*idimp1
  192. xp1 = xcoor(ipv1+1)
  193. yp1 = xcoor(ipv1+2)
  194. zp1 = xcoor(ipv1+3)
  195. ipv2 = (ip2-1)*idimp1
  196. xp2 = xcoor(ipv2+1)
  197. yp2 = xcoor(ipv2+2)
  198. zp2 = xcoor(ipv2+3)
  199. ipv3 = (ip3-1)*idimp1
  200. xp3 = xcoor(ipv3+1)
  201. yp3 = xcoor(ipv3+2)
  202. zp3 = xcoor(ipv3+3)
  203. *
  204. * normale au plan (123)
  205. *
  206. x12 = xp2 - xp1
  207. y12 = yp2 - yp1
  208. z12 = zp2 - zp1
  209. x23 = xp3 - xp2
  210. y23 = yp3 - yp2
  211. z23 = zp3 - zp2
  212. xn = y12*z23 - z12*y23
  213. yn = z12*x23 - x12*z23
  214. zn = x12*y23 - y12*x23
  215. sn = xn*xn + yn*yn + zn*zn
  216. sn = max(sqrt(abs(sn)),xpetit*1d10)
  217. xn = xn/sn
  218. yn = yn/sn
  219. zn = zn/sn
  220. xms(ip1)=xms(ip1)+xn
  221. yms(ip1)=yms(ip1)+yn
  222. zms(ip1)=zms(ip1)+zn
  223. xms(ip2)=xms(ip2)+xn
  224. yms(ip2)=yms(ip2)+yn
  225. zms(ip2)=zms(ip2)+zn
  226. xms(ip3)=xms(ip3)+xn
  227. yms(ip3)=yms(ip3)+yn
  228. zms(ip3)=zms(ip3)+zn
  229. 820 continue
  230. do 821 ip = 1,nbpts
  231. sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)+zms(ip)*zms(ip)
  232. sn = max(sqrt(abs(sn)),xpetit*1d10)
  233. xms(ip)=xms(ip)/sn
  234. yms(ip)=yms(ip)/sn
  235. zms(ip)=zms(ip)/sn
  236. 821 continue
  237. *
  238. * Nombre d'iterations pour la detection de la direction du contact
  239. nbiter=16
  240. *
  241. * Relation du type 3 : noeud-triangle
  242. *-------------------------------------
  243. *
  244. * creation du meleme associe a la relation
  245. *
  246. nbsous=0
  247. nbref =0
  248. nbnn =5
  249. nbelem=incnel
  250. segini meleme
  251. itypel=22
  252. irigel(1,nrigel)=meleme
  253. *
  254. * creation du descriptif commun a toutes les raideurs
  255. *
  256. nligrp = 13
  257. nligrd = nligrp
  258. segini,descr
  259. lisinc(1)='LX '
  260. lisdua(1)='FLX '
  261. noelep(1)=1
  262. noeled(1)=1
  263. do 10 i=2, nligrp, 3
  264. lisinc(i) =modepl(1)
  265. lisinc(i+1)=modepl(2)
  266. lisinc(i+2)=modepl(3)
  267. lisdua(i) =moforc(1)
  268. lisdua(i+1)=moforc(2)
  269. lisdua(i+2)=moforc(3)
  270. noelep(i) =(i+4)/3
  271. noelep(i+1)=noelep(i)
  272. noelep(i+2)=noelep(i)
  273. noeled(i) =noelep(i)
  274. noeled(i+1)=noelep(i)
  275. noeled(i+2)=noelep(i)
  276. 10 continue
  277. segdes,descr
  278. irigel(3,nrigel) = descr
  279. *
  280. * creation du segment xmatri
  281. *
  282. nelrig = incnel
  283. segini,xmatri
  284. irigel(4,nrigel) = xmatri
  285. *
  286. * ce qu'on cree est unilateral
  287. *
  288. irigel(6,nrigel)=1
  289. *
  290. * ce qu'on cree est symetrique
  291. *
  292. irigel(7,nrigel)=0
  293. *
  294. * Nombre d'elements dans la rigidite du type 3
  295. nelri3 = 0
  296.  
  297. * Boucle sur le maillage support du contact(frottement)
  298. *
  299. ip1o=0
  300. ip2o=0
  301. ip3o=0
  302. DO 20 iel=1,nbel1
  303. *
  304. xjr = 0d0
  305. ip1=ipt1.num(2,iel)
  306. ip2=ipt1.num(3,iel)
  307. ip3=ipt1.num(4,iel)
  308. if(ip1.eq.ip1o.and.ip2.eq.ip2o.and.ip3.eq.ip3o) goto 25
  309. * Recuperation des coordonees des noeuds du triangle
  310. ipv = (ip1-1)*idimp1
  311. xp1 = xcoor(ipv+1)
  312. yp1 = xcoor(ipv+2)
  313. zp1 = xcoor(ipv+3)
  314. ipv = (ip2-1)*idimp1
  315. xp2 = xcoor(ipv+1)
  316. yp2 = xcoor(ipv+2)
  317. zp2 = xcoor(ipv+3)
  318. ipv = (ip3-1)*idimp1
  319. xp3 = xcoor(ipv+1)
  320. yp3 = xcoor(ipv+2)
  321. zp3 = xcoor(ipv+3)
  322. * Centre de gravite du triangle
  323. xg = (xp1+xp2+xp3) /3.d0
  324. yg = (yp1+yp2+yp3) /3.d0
  325. zg = (zp1+zp2+zp3) /3.d0
  326. * critere de distance
  327. d1 = (xg-xp1)**2 + (yg-yp1)**2 + (zg-zp1)**2
  328. d2 = (xg-xp2)**2 + (yg-yp2)**2 + (zg-zp2)**2
  329. d3 = (xg-xp3)**2 + (yg-yp3)**2 + (zg-zp3)**2
  330. disto = max(d1,d2,d3)
  331. sqdist = sqrt(disto)
  332. * Triangle un peu plus grand pour les tests
  333. scale=1.00 + xszpre
  334. xp1e=xg+(xp1-xg)*scale
  335. yp1e=yg+(yp1-yg)*scale
  336. zp1e=zg+(zp1-zg)*scale
  337. xp2e=xg+(xp2-xg)*scale
  338. yp2e=yg+(yp2-yg)*scale
  339. zp2e=zg+(zp2-zg)*scale
  340. xp3e=xg+(xp3-xg)*scale
  341. yp3e=yg+(yp3-yg)*scale
  342. zp3e=zg+(zp3-zg)*scale
  343. 25 continue
  344. dist=disto
  345. ip1o=ip1
  346. ip2o=ip2
  347. ip3o=ip3
  348. jp =ipt1.num(5,iel)
  349. ipv = (jp-1)*idimp1
  350. xp = xcoor(ipv+1)
  351. yp = xcoor(ipv+2)
  352. zp = xcoor(ipv+3)
  353. * Verification que pas relation du point sur lui meme
  354. if (jp.eq.ip1) goto 20
  355. if (jp.eq.ip2) goto 20
  356. if (jp.eq.ip3) goto 20
  357. if (MELVA1.ne.0) then
  358. nel1 = min (iel,nelj)
  359. xjr = melva1.velche(nptelj,nel1)
  360. endif
  361. * verification rapide en norme L1 d'abord
  362. if(abs(xp-xg).gt.sqdist*3.) goto 20
  363. if(abs(yp-yg).gt.sqdist*3.) goto 20
  364. if(abs(zp-zg).gt.sqdist*3.) goto 20
  365. *
  366. * Verification si autre point dans la zone de selection
  367. * verif distance au centre de gravite
  368. dp = ((xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2)
  369. if (dp.gt.dist*9) goto 20
  370. C*DBG write(ioimp,*) 'contact test distance ok',dp,d1,d2,d3
  371. dist = SQRT(dist)
  372.  
  373. * verif position par rapport aux cotes
  374.  
  375. * cote 1 2
  376. x12 = xp2 - xp1
  377. y12 = yp2 - yp1
  378. z12 = zp2 - zp1
  379. sv = sqrt(x12**2+y12**2+z12**2)
  380. xv=x12/sv
  381. yv=y12/sv
  382. zv=z12/sv
  383. * normale locale (1 et 2)
  384. xnl=xms(ip1)+xms(ip2)
  385. ynl=yms(ip1)+yms(ip2)
  386. znl=zms(ip1)+zms(ip2)
  387.  
  388. * vecteur reference
  389. xn=y12*znl-z12*ynl
  390. yn=z12*xnl-x12*znl
  391. zn=x12*ynl-y12*xnl
  392. dn = sqrt(xn*xn+yn*yn+zn*zn)
  393. scal = (xp-xp1e)*xv + (yp-yp1e)*yv + (zp-zp1e)*zv
  394. xm = xp1e + scal*xv
  395. ym = yp1e + scal*yv
  396. zm = zp1e + scal*zv
  397. dpm = sqrt(abs((xp-xm)**2+(yp-ym)**2+(zp-zm)**2))
  398. scal = (xp-xm)*xn + (yp-ym)*yn + (zp-zm)*zn
  399. ** if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 1 2',dpm
  400. if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) goto 20
  401. dpm=max(xpetit,dpm)
  402. * write(ioimp,*) 'contact test position 1 ok',
  403. * & scal/(dn*dpm),scal,dn,dpm
  404. *
  405. * cote 2 3
  406. x23 = xp3 - xp2
  407. y23 = yp3 - yp2
  408. z23 = zp3 - zp2
  409. sv = sqrt(x23**2+y23**2+z23**2)
  410. xv=x23/sv
  411. yv=y23/sv
  412. zv=z23/sv
  413. * normale locale (2 et 3)
  414. xnl=xms(ip2)+xms(ip3)
  415. ynl=yms(ip2)+yms(ip3)
  416. znl=zms(ip2)+zms(ip3)
  417.  
  418. * vecteur reference
  419. xn=y23*znl-z23*ynl
  420. yn=z23*xnl-x23*znl
  421. zn=x23*ynl-y23*xnl
  422. dn = sqrt(xn*xn+yn*yn+zn*zn)
  423. scal = (xp-xp2e)*xv + (yp-yp2e)*yv + (zp-zp2e)*zv
  424. xm = xp2e + scal*xv
  425. ym = yp2e + scal*yv
  426. zm = zp2e + scal*zv
  427. dpm = sqrt(abs((xp-xm)**2+(yp-ym)**2+(zp-zm)**2))
  428. scal = (xp-xm)*xn + (yp-ym)*yn + (zp-zm)*zn
  429. ** if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 2 3',dpm
  430. if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) goto 20
  431. dpm=max(xpetit,dpm)
  432. * write(ioimp,*) 'contact test position 2 ok',
  433. * & scal/(dn*dpm),scal,dn,dpm
  434. *
  435. * cote 3 1
  436. x31 = xp1 - xp3
  437. y31 = yp1 - yp3
  438. z31 = zp1 - zp3
  439. sv = sqrt(x31**2+y31**2+z31**2)
  440. xv=x31/sv
  441. yv=y31/sv
  442. zv=z31/sv
  443. * normale locale (3 et 1)
  444. xnl=xms(ip3)+xms(ip1)
  445. ynl=yms(ip3)+yms(ip1)
  446. znl=zms(ip3)+zms(ip1)
  447. * vecteur reference
  448. xn=y31*znl-z31*ynl
  449. yn=z31*xnl-x31*znl
  450. zn=x31*ynl-y31*xnl
  451. dn = sqrt(xn*xn+yn*yn+zn*zn)
  452. scal = (xp-xp3e)*xv + (yp-yp3e)*yv + (zp-zp3e)*zv
  453. xm = xp3e + scal*xv
  454. ym = yp3e + scal*yv
  455. zm = zp3e + scal*zv
  456. dpm = sqrt(abs((xp-xm)**2+(yp-ym)**2+(zp-zm)**2))
  457. scal = (xp-xm)*xn + (yp-ym)*yn + (zp-zm)*zn
  458. ** if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 3 1',dpm
  459. if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) goto 20
  460. dpm=max(xpetit,dpm)
  461. * write(ioimp,*) 'contact test position 2 ok',
  462. * & scal/(dn*dpm),scal,dn,dpm
  463. *
  464. * on a un point ou imposer la relation
  465. C*DBG write(ioimp,*) 'contact potentiel '
  466. * Quelle est la relation ???
  467. *
  468. * direction de la relation = normale au plan (123)
  469. *
  470. * normale reelle
  471. xnr = y12*z23 - z12*y23
  472. ynr = z12*x23 - x12*z23
  473. znr = x12*y23 - y12*x23
  474. dnr = sqrt(xnr*xnr+ynr*ynr+znr*znr)
  475. if (dnr.lt.dist*xszpre.AND.iimpi.ne.0) write(ioimp,*) ' pb dnr'
  476. xnr = xnr / dnr
  477. ynr = ynr / dnr
  478. znr = znr / dnr
  479. * normale ponderee
  480. xn=xms(ip1)+xms(ip2)+xms(ip3)
  481. yn=yms(ip1)+yms(ip2)+yms(ip3)
  482. zn=zms(ip1)+zms(ip2)+zms(ip3)
  483. dn = sqrt(xn*xn+yn*yn+zn*zn)
  484. if (.not.(dn.lt.1d-10).and. .not.(dn.gt.1d-10).AND.iimpi.ne.0)
  485. & write(ioimp,*) 'Prob 4.1 - impo32'
  486. if (abs(dn).le.xpetit.AND.iimpi.ne.0)
  487. & write(ioimp,*) 'Prob 4.2 - impo32'
  488. xn = xn / dn
  489. yn = yn / dn
  490. zn = zn / dn
  491. *
  492. * Distance (jeu) du point jp au plan de contact
  493. * Puis calcul de la projection sur le plan de contact
  494. *
  495. ** write(6,*) 'ip1 xms ',ip1,xms(ip1),yms(ip1),zms(ip1)
  496. ** write(6,*) 'ip2 xms ',ip2,xms(ip2),yms(ip2),zms(ip2)
  497. ** write(6,*) 'ip3 xms ',ip3,xms(ip3),yms(ip3),zms(ip3)
  498. ** write(6,*) 'xn ',xn,yn,zn
  499. iter = 0
  500. 21 continue
  501. xjeuv = (xp-xg)*xnr + (yp-yg)*ynr + (zp-zg)*znr
  502. angn = xn * xnr + yn*ynr + zn*znr
  503. if (angn.lt.0.d0.AND.iimpi.ne.0)
  504. & write(ioimp,*) 'angn negatif ',angn
  505. xjeuv = xjeuv / angn
  506. *
  507. * Recherche de ses coordonnées barycentriques
  508. * qui sont les surfaces des sous triangles
  509. *
  510. xq = xp - xjeuv*xn
  511. yq = yp - xjeuv*yn
  512. zq = zp - xjeuv*zn
  513. xb1 = (yq-yp2)*(zq-zp3)-(zq-zp2)*(yq-yp3)
  514. yb1 = (zq-zp2)*(xq-xp3)-(xq-xp2)*(zq-zp3)
  515. zb1 = (xq-xp2)*(yq-yp3)-(yq-yp2)*(xq-xp3)
  516. xb2 = (yq-yp3)*(zq-zp1)-(zq-zp3)*(yq-yp1)
  517. yb2 = (zq-zp3)*(xq-xp1)-(xq-xp3)*(zq-zp1)
  518. zb2 = (xq-xp3)*(yq-yp1)-(yq-yp3)*(xq-xp1)
  519. xb3 = (yq-yp1)*(zq-zp2)-(zq-zp1)*(yq-yp2)
  520. yb3 = (zq-zp1)*(xq-xp2)-(xq-xp1)*(zq-zp2)
  521. zb3 = (xq-xp1)*(yq-yp2)-(yq-yp1)*(xq-xp2)
  522. b1 = xb1*xn + yb1*yn + zb1*zn
  523. b2 = xb2*xn + yb2*yn + zb2*zn
  524. b3 = xb3*xn + yb3*yn + zb3*zn
  525. bt = b1 + b2 + b3
  526. * normalement pas utile
  527. if (bt.lt.0.d0.AND.iimpi.ne.0) then
  528. write (ioimp,*) ' bt negatif dans impo32 '
  529. print *,'xp1 yp1 zp1',xp1,yp1,zp1
  530. print *,'xp2 yp2 zp2',xp2,yp2,zp2
  531. print *,'xp3 yp3 zp3',xp3,yp3,zp3
  532. print *,'xp yp zp',xp,yp,zp
  533. print *,'xn yn zn',xn,yn,zn
  534. print *,'b1 b2 b3 bt',b1,b2,b3,bt
  535. * goto 20
  536. endif
  537. bsurf = bt / 2
  538. if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10).AND.iimpi.ne.0)
  539. & write(ioimp,*) 'Prob 5.1 - impo32'
  540. if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 5.2 - impo32'
  541. b1 = b1 / bt
  542. b2 = b2 / bt
  543. b3 = b3 / bt
  544. * recalcul de la direction et du jeu en fonction des normales aux sommets
  545. * on arete l'interpolation de la normale au bord de l'element
  546. bb1=max(b1,0.d0)
  547. bb2=max(b2,0.d0)
  548. bb3=max(b3,0.d0)
  549. bm = bb1 + bb2 + bb3
  550. bb1 = bb1 / bm
  551. bb2 = bb2 / bm
  552. bb3 = bb3 / bm
  553. xn = bb1*xms(ip1)+ bb2*xms(ip2)+ bb3*xms(ip3)
  554. yn = bb1*yms(ip1)+ bb2*yms(ip2)+ bb3*yms(ip3)
  555. zn = bb1*zms(ip1)+ bb2*zms(ip2)+ bb3*zms(ip3)
  556. sn = sqrt (xn*xn + yn*yn + zn*zn)
  557. xn=xn/sn
  558. yn=yn/sn
  559. zn=zn/sn
  560. * write (6,*) ' iter b* *n ',iter,b1,b2,b3,xn,yn,zn
  561. * recalcul en fonction de la nouvelle normale
  562. iter=iter+1
  563. if (iter.le.nbiter) goto 21
  564. ** write (6,*) ' iter b* *n ',iter,b1,b2,b3,xn,yn,zn
  565.  
  566.  
  567. * si on a deja traverse, les trois coordonnees barycentriques doivent etre positives
  568. ** if (xjeuv-xjr.lt.-xszpre*dist) then
  569. ** if (b1.lt.-xszpre) goto 20
  570. ** if (b2.lt.-xszpre) goto 20
  571. ** if (b3.lt.-xszpre) goto 20
  572. ** endif
  573. * Si on est en dehors, on projette sur l'arete (ou pas)
  574. ** if (b1.lt.0.d0.or.b2.lt.0.d0.or.b3.lt.0.d0) then
  575. ** xq=xp1*bb1+xp2*bb2+xp3*bb3
  576. ** yq=yp1*bb1+yp2*bb2+yp3*bb3
  577. ** zq=zp1*bb1+zp2*bb2+zp3*bb3
  578. ** xnn=xp-xq
  579. ** ynn=yp-yq
  580. ** znn=zp-zq
  581. ** snn=sqrt(abs(xnn*xnn+ynn*ynn+znn*znn))
  582. ** if (snn.lt.xszpre*dist) write (6,*) ' snn petit ',snn
  583. * sinon on prend la direction reelle
  584. ** if (snn.gt.xszpre*dist) then
  585. ** xn=xnn/sn
  586. ** yn=ynn/sn
  587. ** zn=znn/sn
  588. ** endif
  589. ** endif
  590. xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn
  591. xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn
  592. xjeu3 = (xp-xp3)*xn + (yp-yp3)*yn + (zp-zp3)*zn
  593. xjeuv = xjeu1 * bb1 + xjeu2 * bb2 + xjeu3 * bb3
  594. ** write (6,*) 'b1 b2 b3',b1,b2,b3
  595. ** write (6,*) 'xjeu xjeuv ',xjeu,xjeuv,xjeu1,xjeu2,xjeu3
  596. xjeu = xjeuv - xjr
  597. ** write (6,*) ' normale ',xn,yn,zn,' jeu ',xjeu1,xjeu2,xjeu3
  598. * verif bon cote (a la tolerance pres)
  599. ** write (6,*) ' xjeu ',xjeu,dist,iel
  600. if (xjeu.lt.-3*dist) goto 20
  601. if (xjeu.gt.3*dist) goto 20
  602. * verif compatibilite avec la normale au poin impactant
  603. if (xms(jp)*xn+yms(jp)*yn+zms(jp)*zn.gt. -0.0d0) then
  604. ** write (6,*) ' impo32 normales incompatibes 1 ',
  605. ** > jp,xjeu,dpm,dist
  606. goto 20
  607. endif
  608.  
  609.  
  610.  
  611. C*DBG write(ioimp,*) ' b1 b2 b3 ',b1,b2,b3
  612. * creation du support du multiplicateur de lagrange si necessaire
  613. ** nbpts = xcoor(/1)/(idim+1)
  614. if (ipt1.num(1,iel).eq.ilgni) then
  615. nbpts = nbpts + 1
  616. segadj mcoord
  617. segact ipt1*mod
  618. ipt1.num(1,iel)=nbpts
  619. if (ipt6.ne.0) ipt6.num(1,iel)=nbpts
  620. endif
  621. segact mcoord*mod
  622. xcoor(nbpts*(idim+1)) = dist
  623. 1954 continue
  624. * points a l'exterieur?
  625. * on met une ponderation et un rayon d'acceptation
  626. pond = (1D4+1) - 1D4 * bm
  627. if (pond.le.0.d0) goto 20
  628. pond=max(pond,0.d0)
  629. pond=min(pond,1.d0)
  630.  
  631. ** if (xjeu.lt.0.d0)
  632. ** > write (6,*) 'pt traverse',jp,b1,b2,b3,xjeuv-xjr,dist
  633.  
  634. *
  635. * Ajout d'une relation noeud-triangle
  636. *
  637. nelri3 = nelri3 + 1
  638. nelch = nelch + 1
  639. *
  640. * on ajuste les differents segments si necesaire
  641. *
  642. if (nelri3.gt.nelrig) then
  643. nelrig = nelrig + incnel
  644. segadj,xmatri
  645. nbelem = nbelem + incnel
  646. nbnn = 5
  647. segadj,meleme
  648. nbnn = 1
  649. segadj,ipt7
  650. n = n + incnel
  651. segadj,mpoval
  652. endif
  653. *
  654. * Mise a jour du meleme
  655. num(1,nelri3) = ipt1.num(1,iel)
  656. num(2,nelri3) = ip1
  657. num(3,nelri3) = ip2
  658. num(4,nelri3) = ip3
  659. num(5,nelri3) = jp
  660. *
  661. * Mise a jour de xmatri
  662. * lambda
  663. re( 1,1,nelri3) = 0d0
  664. * ip1
  665. re( 2,1,nelri3) = xn * bb1 * bsurf * pond
  666. re( 3,1,nelri3) = yn * bb1 * bsurf * pond
  667. re( 4,1,nelri3) = zn * bb1 * bsurf * pond
  668. * ip2
  669. re( 5,1,nelri3) = xn * bb2 * bsurf * pond
  670. re( 6,1,nelri3) = yn * bb2 * bsurf * pond
  671. re( 7,1,nelri3) = zn * bb2 * bsurf * pond
  672. * ip3
  673. re( 8,1,nelri3) = xn * bb3 * bsurf * pond
  674. re( 9,1,nelri3) = yn * bb3 * bsurf * pond
  675. re(10,1,nelri3) = zn * bb3 * bsurf * pond
  676. * jp
  677. re(11,1,nelri3) = -xn * bsurf * pond
  678. re(12,1,nelri3) = -yn * bsurf * pond
  679. re(13,1,nelri3) = -zn * bsurf * pond
  680. * on transpose
  681. do 30 ic = 2, nligrp
  682. re(1,ic,nelri3) = re(ic,1,nelri3)
  683. 30 continue
  684. * le reste est nul
  685. *
  686. * remplissage du champoint de depi (jeu)
  687. *
  688. ipt7.num(1,nelch) = ipt1.num(1,iel)
  689. vpocha(nelch,1) = xjeu * bsurf * pond
  690. IF (MELVA2.NE.0) THEN
  691. NEL2 = min (iel,NELA)
  692. VPOCHA(nelch,2) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*
  693. > bsurf * pond
  694. ENDIF
  695. * write(ioimp,*) ' jeu ',xjeu,ip1,ip2,ip3,jp,b1,b2,b3,xn,yn,zn
  696. *
  697. 20 CONTINUE
  698. *
  699. * Ajustement au plus juste puis desactivation des segments lies
  700. * a la rigidite du type 3
  701. if (nelri3.ne.nelrig) then
  702. nelrig = nelri3
  703. segadj,xmatri
  704. nbelem = nelri3
  705. nbnn = 5
  706. segadj,meleme
  707. endif
  708. segdes,xmatri
  709. *
  710. *pv write(ioimp,*) ' nb relation type 3 ',nelri3
  711. ** write(6,fmt='(10i8)') (num(5,nelr),nelr=1,nelri3)
  712.  
  713.  
  714.  
  715.  
  716.  
  717. C*DBG if (nelri3.eq.0) irigel(6,nrigel)=0
  718. *
  719. * Fin du traitement de la formulation standard/forte
  720. * ----------------------------------------------------
  721. 300 CONTINUE
  722. * Destruction des segments locaux devenus inutiles
  723. segsup mfopa3
  724.  
  725. GOTO 1000
  726.  
  727. *=======================================================================
  728. *= Formulation "faible" 3D : relation triangle-triangle
  729. *=======================================================================
  730. * Element de contact 3D a 7 noeuds (1+6)
  731. 500 CONTINUE
  732. *
  733. * Petit segment de travail
  734. segini,mfaible
  735. *
  736. * Relation du type 0 : triangle-triangle (faible)
  737. *---------------------------------------
  738. *
  739. * Creation du meleme associe a la relation
  740. *
  741. nbnn = 7
  742. nbelem = incnel
  743. nbsous = 0
  744. nbref = 0
  745. segini,meleme
  746. itypel = 22
  747. irigel(1,nrigel) = meleme
  748. *
  749. * Creation du descriptif commun a toutes les raideurs
  750. *
  751. nligrp = 19
  752. nligrd = nligrp
  753. segini,descr
  754. lisinc(1) = 'LX '
  755. lisdua(1) = 'FLX '
  756. noelep(1) = 1
  757. noeled(1) = 1
  758. do 510 i = 2, nligrp, 3
  759. lisinc(i ) = modepl(1)
  760. lisinc(i+1) = modepl(2)
  761. lisinc(i+2) = modepl(3)
  762. lisdua(i ) = moforc(1)
  763. lisdua(i+1) = moforc(2)
  764. lisdua(i+2) = moforc(3)
  765. noelep(i ) = (i+5)/3
  766. noelep(i+1) = noelep(i)
  767. noelep(i+2) = noelep(i)
  768. noeled(i ) = noelep(i)
  769. noeled(i+1) = noelep(i)
  770. noeled(i+2) = noelep(i)
  771. 510 continue
  772. segdes,descr
  773. irigel(3,nrigel) = descr
  774. *
  775. * creation du segment xmatri
  776. *
  777. nelrig = incnel
  778. segini,xmatri
  779. irigel(4,nrigel) = xmatri
  780. *
  781. * ce qu'on cree est unilateral
  782. *
  783. irigel(6,nrigel) = 1
  784. *
  785. * ce qu'on cree est symetrique
  786. *
  787. irigel(7,nrigel) = 0
  788. *
  789. * Nombre d'elements crees dans meleme=irigel(1,nrigel), ipt7 et mpoval
  790. nelri0 = 0
  791. *
  792. * Boucle sur les elements du maillage de contact/frottement "faible"
  793. *
  794. DO 520 iel = 1, nbel1
  795. *
  796. xjr = 0d0
  797. ip1 = ipt1.num(2,iel)
  798. ip2 = ipt1.num(3,iel)
  799. ip3 = ipt1.num(4,iel)
  800. jp1 = ipt1.num(5,iel)
  801. jp2 = ipt1.num(6,iel)
  802. jp3 = ipt1.num(7,iel)
  803. if (MELVA1.ne.0) then
  804. nel1 = min (iel,nelj)
  805. xjr = melva1.velche(nptelj,nel1)
  806. endif
  807. C*DBG write(ioimp,*) iel,ip1,ip2,ip3,jp1,jp2,jp3,ipt1.num(1,iel)
  808. *
  809. * Verification (provisoire) qu'il n'y a pas de noeud double
  810. if (ip1.eq.jp1) goto 520
  811. if (ip1.eq.jp2) goto 520
  812. if (ip1.eq.jp3) goto 520
  813. if (ip2.eq.jp1) goto 520
  814. if (ip2.eq.jp2) goto 520
  815. if (ip2.eq.jp3) goto 520
  816. if (ip3.eq.jp1) goto 520
  817. if (ip3.eq.jp2) goto 520
  818. if (ip3.eq.jp3) goto 520
  819. *
  820. * Definition du triangle T1(ip1,ip2,ip3)
  821. * - Coordonnees des noeuds de T1
  822. ipv = (ip1-1)*idimp1
  823. xp1 = xcoor(ipv+1)
  824. yp1 = xcoor(ipv+2)
  825. zp1 = xcoor(ipv+3)
  826. ipv = (ip2-1)*idimp1
  827. xp2 = xcoor(ipv+1)
  828. yp2 = xcoor(ipv+2)
  829. zp2 = xcoor(ipv+3)
  830. ipv = (ip3-1)*idimp1
  831. xp3 = xcoor(ipv+1)
  832. yp3 = xcoor(ipv+2)
  833. zp3 = xcoor(ipv+3)
  834. * - Barycentre de T1
  835. xgT1 = (xp1 + xp2 + xp3) * X1s3
  836. ygT1 = (yp1 + yp2 + yp3) * X1s3
  837. zgT1 = (zp1 + zp2 + zp3) * X1s3
  838. * - Normale au triangle T1
  839. xnT1 = (yp1-yp2)*(zp2-zp3) - (zp1-zp2)*(yp2-yp3)
  840. ynT1 = (zp1-zp2)*(xp2-xp3) - (xp1-xp2)*(zp2-zp3)
  841. znT1 = (xp1-xp2)*(yp2-yp3) - (yp1-yp2)*(xp2-xp3)
  842. dnT1 = SQRT(xnT1*xnT1+ynT1*ynT1+znT1*znT1)
  843. IF (.not.(dnT1.lt.1d-10) .and. .not.(dnT1.gt.1d-10).AND.
  844. & iimpi.ne.0) write(ioimp,*) 'FAIBle - 1.1 - impo32'
  845. IF (abs(dnT1).le.xpetit) write(ioimp,*) 'FAIBle - 1.2 - impo32'
  846. xnT1 = xnT1 / dnT1
  847. ynT1 = ynT1 / dnT1
  848. znT1 = znT1 / dnT1
  849. C*DBG write(ioimp,*) 'Normale au plan T1',xnT1,ynT1,znT1,dnT1
  850. * - Rayon de la sphere "englobante" de l'element (centre=barycentre)
  851. d1 = (xgT1-xp1)**2 + (ygT1-yp1)**2 + (zgT1-zp1)**2
  852. d2 = (xgT1-xp2)**2 + (ygT1-yp2)**2 + (zgT1-zp2)**2
  853. d3 = (xgT1-xp3)**2 + (ygT1-yp3)**2 + (zgT1-zp3)**2
  854. RayT1 = SQRT(max(d1,d2,d3))
  855. *
  856. * Definition du triangle T2(jp1,jp2,jp3)
  857. * - Coordonnees des noeuds de T2
  858. ipv = (jp1-1)*idimp1
  859. xq1 = xcoor(ipv+1)
  860. yq1 = xcoor(ipv+2)
  861. zq1 = xcoor(ipv+3)
  862. ipv = (jp2-1)*idimp1
  863. xq2 = xcoor(ipv+1)
  864. yq2 = xcoor(ipv+2)
  865. zq2 = xcoor(ipv+3)
  866. ipv = (jp3-1)*idimp1
  867. xq3 = xcoor(ipv+1)
  868. yq3 = xcoor(ipv+2)
  869. zq3 = xcoor(ipv+3)
  870. * - Barycentre de T2
  871. xgT2 = (xq1 + xq2 + xq3) * X1s3
  872. ygT2 = (yq1 + yq2 + yq3) * X1s3
  873. zgT2 = (zq1 + zq2 + zq3) * X1s3
  874. * - Normale au triangle T2
  875. xnT2 = (yq1-yq2)*(zq2-zq3) - (zq1-zq2)*(yq2-yq3)
  876. ynT2 = (zq1-zq2)*(xq2-xq3) - (xq1-xq2)*(zq2-zq3)
  877. znT2 = (xq1-xq2)*(yq2-yq3) - (yq1-yq2)*(xq2-xq3)
  878. dnT2 = SQRT(xnT2*xnT2+ynT2*ynT2+znT2*znT2)
  879. IF (.not.(dnT2.lt.1d-10) .and. .not.(dnT2.gt.1d-10).AND.
  880. & iimpi.ne.0) write(ioimp,*) 'FAIBle - 2.1 - impo32'
  881. IF (abs(dnT2).le.xpetit.AND.iimpi.ne.0)
  882. & write(ioimp,*) 'FAIBle - 2.2 - impo32'
  883. xnT2 = xnT2 / dnT2
  884. ynT2 = ynT2 / dnT2
  885. znT2 = znT2 / dnT2
  886. C*DBG write(ioimp,*) 'Normale au plan T2',xnT2,ynT2,znT2,dnT2
  887. * - Rayon de la sphere "englobante" de l'element (centre=barycentre)
  888. d1 = (xgT2-xq1)**2 + (ygT2-yq1)**2 + (zgT2-zq1)**2
  889. d2 = (xgT2-xq2)**2 + (ygT2-yq2)**2 + (zgT2-zq2)**2
  890. d3 = (xgT2-xq3)**2 + (ygT2-yq3)**2 + (zgT2-zq3)**2
  891. RayT2 = SQRT(max(d1,d2,d3))
  892.  
  893. * Orientation respective correcte des 2 normales
  894. * Les triangles T1 et T2 etant decrits dans le meme sens ("prisme"),
  895. * le produit scalaire de leur normale doit etre positif ou nul !
  896. scal = xnT1 * xnT2 + ynT1 * ynT2 + znT1 * znT2
  897. C*DBG write(ioimp,*) 'Prod.scal des normales',scal
  898. if (scal.lt.0.) goto 520
  899. *
  900. * Distance entre les centres de gravite
  901. * Les spheres "englobantes" doivent s'intersecter pour contact potentiel
  902. dist = SQRT((xgT2-xgT1)**2 + (ygT2-ygT1)**2 + (zgT2-zgT1)**2)
  903. IF (dist*4.GT.RayT1+RayT2) goto 520
  904. *
  905. * Critere d'acceptation de l'element de contact
  906. *
  907. * Definition du plan de contact (surface intermediaire)
  908. * - Normale (unitaire) du plan de contact
  909. xnC = xnT1 + xnT2
  910. ynC = ynT1 + ynT2
  911. znC = znT1 + znT2
  912. dnC = SQRT(xnC*xnC + ynC*ynC + znC*znC)
  913. IF (.not.(dnC.lt.1d-10).and. .not.(dnC.GT.1D-10).AND.iimpi.ne.0)
  914. & write(ioimp,*) 'FAIBle - 3.1 - impo32'
  915. IF (abs(dnC).le.xpetit.AND.iimpi.ne.0)
  916. & write(ioimp,*) 'FAIBle - 3.2 - impo32'
  917. xnC = xnC / dnC
  918. ynC = ynC / dnC
  919. znC = znC / dnC
  920. C*DBG write(ioimp,*) 'Normale au plan Contact',xnC,ynC,znC,dnC
  921.  
  922. * - Point definissant le plan de contact
  923. * = milieu des barycentres des triangles T1 et T2
  924. xgC = (xgT1 + xgT2) * 0.5
  925. ygC = (ygT1 + ygT2) * 0.5
  926. zgC = (zgT1 + zgT2) * 0.5
  927. * - "Distance" de reference au plan de contact
  928. dgC = xgC*xnC + ygC*ynC + zgC*znC
  929. *
  930. * Distance "signee" des triangles T1 et T2 au plan de contact
  931. dp1C = (xp1*xnC + yp1*ynC + zp1*znC) - dgC
  932. dp2C = (xp2*xnC + yp2*ynC + zp2*znC) - dgC
  933. dp3C = (xp3*xnC + yp3*ynC + zp3*znC) - dgC
  934. dq1C = (xq1*xnC + yq1*ynC + zq1*znC) - dgC
  935. dq2C = (xq2*xnC + yq2*ynC + zq2*znC) - dgC
  936. dq3C = (xq3*xnC + yq3*ynC + zq3*znC) - dgC
  937. * Projection des triangles T1 et T2 sur le plan de contact
  938. * T1C = triangle T1 projete sur plan de Contact
  939. cPT1C(1,1) = xp1 - dp1C * xnC
  940. cPT1C(1,2) = yp1 - dp1C * ynC
  941. cPT1C(1,3) = zp1 - dp1C * znC
  942. cPT1C(2,1) = xp2 - dp2C * xnC
  943. cPT1C(2,2) = yp2 - dp2C * ynC
  944. cPT1C(2,3) = zp2 - dp2C * znC
  945. cPT1C(3,1) = xp3 - dp3C * xnC
  946. cPT1C(3,2) = yp3 - dp3C * ynC
  947. cPT1C(3,3) = zp3 - dp3C * znC
  948. * T2C = triangle T2 projete sur plan de Contact
  949. cPT2C(1,1) = xq1 - dq1C * xnC
  950. cPT2C(1,2) = yq1 - dq1C * ynC
  951. cPT2C(1,3) = zq1 - dq1C * znC
  952. cPT2C(2,1) = xq2 - dq2C * xnC
  953. cPT2C(2,2) = yq2 - dq2C * ynC
  954. cPT2C(2,3) = zq2 - dq2C * znC
  955. cPT2C(3,1) = xq3 - dq3C * xnC
  956. cPT2C(3,2) = yq3 - dq3C * ynC
  957. cPT2C(3,3) = zq3 - dq3C * znC
  958. *
  959. * On determine quelle est la composante maximale de la normale pour
  960. * projeter les triangles T1C et T2C sur un plan parallele aux axes
  961. * ("le plus orthogonal" a cette normale), ce qui maximise leur aire.
  962. r_x = abs(xnC)
  963. r_y = abs(ynC)
  964. r_z = abs(znC)
  965. if (r_x.gt.r_y) then
  966. inC = 1
  967. if (r_x.lt.r_z) inC = 3
  968. else
  969. inC = 2
  970. if (r_y.lt.r_z) inC = 3
  971. endif
  972. * Projection des triangles T1C et T2C sur ce plan // aux axes (A)
  973. * T1A = triangle T1C projete sur ce plan (A)
  974. * T2A = triangle T2C projete sur ce plan (A)
  975. * On prend soin d'orienter les triangles T1A et T2A dans le sens
  976. * direct (en tenant compte du signe de la composante normale !)
  977. GOTO (531,532,533), inC
  978. call erreur(5)
  979. return
  980. * Normale NC selon "x" : A =(y,z)
  981. 531 if (xnC.lt.0.) inC = -inC
  982. inX = 2
  983. inY = 3
  984. goto 534
  985. * Normale NC selon "y" : A =(x,z)
  986. 532 if (ynC.lt.0.) inC = -inC
  987. inX = 1
  988. inY = 3
  989. goto 534
  990. * Normale NC selon "z" : A =(x,y)
  991. 533 if (znC.lt.0.) inC = -inC
  992. inX = 1
  993. inY = 2
  994. goto 534
  995. 534 continue
  996. cPT1A(1,1) = cPT1C(1,inX)
  997. cPT1A(1,2) = cPT1C(1,inY)
  998. cPT2A(1,1) = cPT2C(1,inX)
  999. cPT2A(1,2) = cPT2C(1,inY)
  1000. if (inC.gt.0) then
  1001. cPT1A(2,1) = cPT1C(2,inX)
  1002. cPT1A(2,2) = cPT1C(2,inY)
  1003. cPT2A(2,1) = cPT2C(2,inX)
  1004. cPT2A(2,2) = cPT2C(2,inY)
  1005. cPT1A(3,1) = cPT1C(3,inX)
  1006. cPT1A(3,2) = cPT1C(3,inY)
  1007. cPT2A(3,1) = cPT2C(3,inX)
  1008. cPT2A(3,2) = cPT2C(3,inY)
  1009. else
  1010. cPT1A(2,1) = cPT1C(3,inX)
  1011. cPT1A(2,2) = cPT1C(3,inY)
  1012. cPT2A(2,1) = cPT2C(3,inX)
  1013. cPT2A(2,2) = cPT2C(3,inY)
  1014. cPT1A(3,1) = cPT1C(2,inX)
  1015. cPT1A(3,2) = cPT1C(2,inY)
  1016. cPT2A(3,1) = cPT2C(2,inX)
  1017. cPT2A(3,2) = cPT2C(2,inY)
  1018. endif
  1019. cPT1A(4,1) = cPT1A(1,1)
  1020. cPT1A(4,2) = cPT1A(1,2)
  1021. * Calcul des cotes des triangles T1A et T2A (12,23,31)
  1022. * Utile pour connaitre ensuite la normale aux cotes des triangles
  1023. vT1A(1,1) = cPT1A(2,1) - cPT1A(1,1)
  1024. vT1A(1,2) = cPT1A(2,2) - cPT1A(1,2)
  1025. vT1A(2,1) = cPT1A(3,1) - cPT1A(2,1)
  1026. vT1A(2,2) = cPT1A(3,2) - cPT1A(2,2)
  1027. vT1A(3,1) = cPT1A(1,1) - cPT1A(3,1)
  1028. vT1A(3,2) = cPT1A(1,2) - cPT1A(3,2)
  1029. C*nu vT2A(1,1) = cPT2A(2,1) - cPT2A(1,1)
  1030. vT2A(1,2) = cPT2A(2,2) - cPT2A(1,2)
  1031. C*nu vT2A(2,1) = cPT2A(3,1) - cPT2A(2,1)
  1032. vT2A(2,2) = cPT2A(3,2) - cPT2A(2,2)
  1033. C*nu vT2A(3,1) = cPT2A(1,1) - cPT2A(3,1)
  1034. vT2A(3,2) = cPT2A(1,2) - cPT2A(3,2)
  1035. * Calcul de la surface de chacun des triangles T1A et T2A
  1036. * (en fait, on calcule le double, mais par la suite on ne considere que
  1037. * des rapports de surfaces de triangle)
  1038. SuT1A = (cPT1A(2,1)+cPT1A(1,1)) * vT1A(1,2)
  1039. & + (cPT1A(3,1)+cPT1A(2,1)) * vT1A(2,2)
  1040. & + (cPT1A(1,1)+cPT1A(3,1)) * vT1A(3,2)
  1041. SuT2A = (cPT2A(2,1)+cPT2A(1,1)) * vT2A(1,2)
  1042. & + (cPT2A(3,1)+cPT2A(2,1)) * vT2A(2,2)
  1043. & + (cPT2A(1,1)+cPT2A(3,1)) * vT2A(3,2)
  1044. C*DBG write(ioimp,*) 'Surfaces T1A - T2A :', 0.5*SuT1A,0.5*SuT2A
  1045. C*DBG if (SuT1A.gt.100.*SuT2A .or. SuT2A.gt.100.*SuT1A)
  1046. C*DBG & write(ioimp,*) 'Rapport des surfaces tres important !'
  1047.  
  1048. * On initialise l'intersection des triangles avec T2A
  1049. nPIn0 = 3
  1050. do 540 i = 1, nPIn0
  1051. cPIn0(i,1) = cPT2A(i,1)
  1052. cPIn0(i,2) = cPT2A(i,2)
  1053. 540 continue
  1054. * Determination de l'intersection des triangles T1A et T2A
  1055. * en regardant progressivement l'intersection avec les cotes de T1A
  1056. * Note : On utilise le fait que les polygones traites sont convexes !
  1057. DO 550 i = 1, 3
  1058. vN1x = -vT1A(i,2)
  1059. vN1y = vT1A(i,1)
  1060. scal = vN1x * cPT1A(i,1) + vN1y * cPT1A(i,2)
  1061. *
  1062. ipos = 0
  1063. ineg = 0
  1064. jpos = 0
  1065. jneg = 0
  1066. do 551 j = 1, nPIn0
  1067. test(j) = (vN1x * cPIn0(j,1) + vN1y * cPIn0(j,2)) - scal
  1068. if (test(j).gt.1.d-10) then
  1069. * if (test(j).gt.0.) then
  1070. ipos = ipos + 1
  1071. if (jneg.ne.0 .and. jpos.eq.0) jpos = j
  1072. else if (test(j).lt.-1.d-10) then
  1073. * else if (test(j).lt.0.) then
  1074. ineg = ineg + 1
  1075. if (jneg.eq.0) jneg = j
  1076. else
  1077. test(j) = 0.
  1078. endif
  1079. 551 continue
  1080. *
  1081. * 1) Cas ou ipos = 0 :
  1082. * 1.1) il n'y pas d'intersection (ineg=nPIn0)
  1083. * 1.2) l'intersection se limite a 1 point (ineg=nPIn0-1)
  1084. * 1.3) l'intersection a 1 segment sur une frontiere (ineg=nPIn0-2)
  1085. * La surface correspondante est nulle et il n'y a donc pas de contact !
  1086. if (ipos.eq.0) then
  1087. C*DBG write(ioimp,*) 'Intersection a ',nPIn0-ineg,' points'
  1088. goto 520
  1089. endif
  1090. * 2) Cas ou ipos > 0 : Il y a une intersection a calculer
  1091. * 2.1) ineg = 0 : Tous les points sont "du bon cote" et l'intersection
  1092. * correspond au polygone de depart !
  1093. if (ineg.eq.0) then
  1094. C*DBG write(ioimp,*) 'Intersection du bon cote - rien a faire'
  1095. goto 550
  1096. endif
  1097. * 2.2) ineg > 0 : Il y a deux intersections a determiner car il y a
  1098. * au moins un point du "mauvais cote"
  1099. jpos = max(1,jpos)
  1100. jpr = jpos - 1
  1101. if (jpos.eq.1) jpr = nPIn0
  1102. nPIn = 1
  1103. if (test(jpr).lt.0.) then
  1104. r_z = test(jpos) / (test(jpos) - test(jpr))
  1105. cPIn(nPIn,1) = cPIn0(jpos,1)
  1106. & + r_z*(cPIn0(jpr,1) - cPIn0(jpos,1))
  1107. cPIn(nPIn,2) = cPIn0(jpos,2)
  1108. & + r_z*(cPIn0(jpr,2) - cPIn0(jpos,2))
  1109. else
  1110. cPIn(nPIn,1) = cPIn0(jpr,1)
  1111. cPIn(nPIn,2) = cPIn0(jpr,2)
  1112. endif
  1113. do 552 j = 1, ipos
  1114. nPIn = nPIn + 1
  1115. cPIn(nPIn,1) = cPIn0(jpos,1)
  1116. cPIn(nPIn,2) = cPIn0(jpos,2)
  1117. jpr = jpos
  1118. jpos = jpos + 1
  1119. if (jpos.gt.nPIn0) jpos = 1
  1120. 552 continue
  1121. nPIn = nPIn + 1
  1122. if (test(jpos).lt.0.) then
  1123. r_z = test(jpr) / (test(jpr) - test(jpos))
  1124. cPIn(nPIn,1) = cPIn0(jpr,1)
  1125. & + r_z*(cPIn0(jpos,1) - cPIn0(jpr,1))
  1126. cPIn(nPIn,2) = cPIn0(jpr,2)
  1127. & + r_z*(cPIn0(jpos,2) - cPIn0(jpr,2))
  1128. else
  1129. cPIn(nPIn,1) = cPIn0(jpos,1)
  1130. cPIn(nPIn,2) = cPIn0(jpos,2)
  1131. endif
  1132. * Mise a jour cPIn0 pour la suite du traitement des intersections...
  1133. nPIn0 = nPIn
  1134. do 554 j = 1, nPIn0
  1135. cPIn0(j,1) = cPIn(j,1)
  1136. cPIn0(j,2) = cPIn(j,2)
  1137. 554 continue
  1138. C*DBG-F write(ioimp,*) 'Intersection ',i,' a ',nPIn0,' points'
  1139. C*DBG-F do j = 1, nPIn0
  1140. C*DBG-F write(ioimp,*) ' ',j,nbpts0+6+j,cPIn0(j,1),cPIn0(j,2)
  1141. C*DBG-F enddo
  1142. *
  1143. 550 CONTINUE
  1144. *
  1145. * Calcul de la surface de l'intersection et de son centre de gravite
  1146. r_z = cPIn0(nPIn0,1)*cPIn0(1,2) - cPIn0(1,1)*cPIn0(nPIn0,2)
  1147. SuPIn = r_z
  1148. xGIn = (cPIn0(nPIn0,1)+cPIn0(1,1)) * r_z
  1149. yGIn = (cPIn0(nPIn0,2)+cPIn0(1,2)) * r_z
  1150. do 560 i = 1, nPIn0-1
  1151. j = i + 1
  1152. r_z = cPIn0(i,1)*cPIn0(j,2) - cPIn0(j,1)*cPIn0(i,2)
  1153. SuPIn = SuPIn + r_z
  1154. xGIn = xGIn + (cPIn0(i,1)+cPIn0(j,1)) * r_z
  1155. yGIn = yGIn + (cPIn0(i,2)+cPIn0(j,2)) * r_z
  1156. 560 continue
  1157. if (SuPIn .lt. 1.E-7*max(SuT1A,SuT2A) ) then
  1158. C*DBG write(ioimp,*) 'Intersection a une surface negligeable !'
  1159. goto 520
  1160. endif
  1161. r_z = X1s3 / SupIn
  1162. SuPIn = 0.5 * SupIn
  1163. xGIn = xGIn * r_z
  1164. yGIn = yGIn * r_z
  1165. C*DBG write(ioimp,*) 'Surface Intersection :',SuPIn
  1166. *
  1167. * Calcul des coordonnees barycentriques du centre de gravite
  1168. * pour le triangle T1A
  1169. xpip1 = xGIn - cPT1A(1,1)
  1170. ypip1 = yGIn - cPT1A(1,2)
  1171. xpip2 = xGIn - cPT1A(2,1)
  1172. ypip2 = yGIn - cPT1A(2,2)
  1173. xpip3 = xGIn - cPT1A(3,1)
  1174. ypip3 = yGIn - cPT1A(3,2)
  1175. b1T1 = xpip2 * ypip3 - ypip2 * xpip3
  1176. b2T1 = xpip3 * ypip1 - ypip3 * xpip1
  1177. b3T1 = xpip1 * ypip2 - ypip1 * xpip2
  1178. bt = b1T1 + b2T1 + b3T1
  1179. if (abs(bt-SuT1A) .gt. 1.E-3)
  1180. & write(ioimp,*) 'Prob bt-SuT1A',bt,SuT1A
  1181. b1T1 = b1T1 / SuT1A
  1182. b2T1 = b2T1 / SuT1A
  1183. b3T1 = b3T1 / SuT1A
  1184. *
  1185. * pour le triangle T2A
  1186. xpip1 = xGIn - cPT2A(1,1)
  1187. ypip1 = yGIn - cPT2A(1,2)
  1188. xpip2 = xGIn - cPT2A(2,1)
  1189. ypip2 = yGIn - cPT2A(2,2)
  1190. xpip3 = xGIn - cPT2A(3,1)
  1191. ypip3 = yGIn - cPT2A(3,2)
  1192. b1T2 = xpip2 * ypip3 - ypip2 * xpip3
  1193. b2T2 = xpip3 * ypip1 - ypip3 * xpip1
  1194. b3T2 = xpip1 * ypip2 - ypip1 * xpip2
  1195. bt = b1T2 + b2T2 + b3T2
  1196. if (abs(bt-SuT2A) .gt. 1.E-3)
  1197. & write(ioimp,*) 'Prob bt-SuT2A',bt,SuT2A
  1198. b1T2 = b1T2 / SuT2A
  1199. b2T2 = b2T2 / SuT2A
  1200. b3T2 = b3T2 / SuT2A
  1201. *
  1202. * Calcul du jeu
  1203. xjeu = (b1T2*dq1C + b2T2*dq2C + b3T2*dq3C)
  1204. & - (b1T1*dp1C + b2T1*dp2C + b3T1*dp3C)
  1205. C*DBG write(ioimp,*) 'Jeu =',xjeu
  1206. *
  1207. * Ajout d'une relation triangle-triangle
  1208. *
  1209. nelri0 = nelri0 + 1
  1210. nelch = nelch + 1
  1211. *
  1212. * on ajuste les differents segments si necesaire
  1213. *
  1214. if (nelri0.gt.nelrig) then
  1215. nelrig = nelrig + incnel
  1216. segadj,xmatri
  1217. nbelem = nbelem + incnel
  1218. nbnn = 7
  1219. segadj,meleme
  1220. nbnn = 1
  1221. segadj,ipt7
  1222. n = n + incnel
  1223. segadj,mpoval
  1224. endif
  1225. *
  1226. * Mise a jour du meleme
  1227. *
  1228. if (ipt1.num(1,iel).eq.ILGNI) then
  1229. segact ipt1*mod
  1230. nbpts=nbpts+1
  1231. segadj mcoord
  1232. ipt1.num(1,iel)=nbpts
  1233. if (ipt6.ne.0) ipt6.num(1,iel)=nbpts
  1234. endif
  1235. num(1,nelri0) = ipt1.num(1,iel)
  1236. num(2,nelri0) = ip1
  1237. num(3,nelri0) = ip2
  1238. num(4,nelri0) = ip3
  1239. num(5,nelri0) = jp1
  1240. num(6,nelri0) = jp2
  1241. num(7,nelri0) = jp3
  1242. icolor(nelri0) = 1
  1243. *$
  1244. * Mise a jour de xmatri
  1245. *
  1246. * lambda
  1247. re( 1,1,nelri0) = 0.d0
  1248. * ip1
  1249. re( 2,1,nelri0) = xnC * b1T1 * SuPIn
  1250. re( 3,1,nelri0) = ynC * b1T1 * SuPIn
  1251. re( 4,1,nelri0) = znC * b1T1 * SuPIn
  1252. * ip2
  1253. re( 5,1,nelri0) = xnC * b2T1 * SuPIn
  1254. re( 6,1,nelri0) = ynC * b2T1 * SuPin
  1255. re( 7,1,nelri0) = znC * b2T1 * SuPIn
  1256. * ip3
  1257. re( 8,1,nelri0) = xnC * b3T1 * SuPIn
  1258. re( 9,1,nelri0) = ynC * b3T1 * SuPIn
  1259. re(10,1,nelri0) = znC * b3T1 * SuPIn
  1260. * jp1
  1261. re(11,1,nelri0) = -xnC * b1T2 * SuPIn
  1262. re(12,1,nelri0) = -ynC * b1T2 * SuPIn
  1263. re(13,1,nelri0) = -znC * b1T2 * SuPIn
  1264. * jp2
  1265. re(14,1,nelri0) = -xnC * b2T2 * SuPIn
  1266. re(15,1,nelri0) = -ynC * b2T2 * SuPIn
  1267. re(16,1,nelri0) = -znC * b2T2 * SuPIn
  1268. * jp3
  1269. re(17,1,nelri0) = -xnC * b3T2 * SuPIn
  1270. re(18,1,nelri0) = -ynC * b3T2 * SuPIn
  1271. re(19,1,nelri0) = -znC * b3T2 * SuPIn
  1272. * on transpose
  1273. do 580 ic = 2, nligrp
  1274. re(1,ic,nelri0) = re(ic,1,nelri0)
  1275. 580 continue
  1276. * le reste est nul
  1277. *
  1278. * remplissage du champoint de depi (jeu)
  1279. *
  1280. ipt7.num(1,nelch) = ipt1.num(1,iel)
  1281. vpocha(nelch,1) = (xjeu - xjr) * SuPIn
  1282. IF (MELVA2.NE.0) THEN
  1283. NEL2 = min (iel,NELA)
  1284. VPOCHA(nelch,2) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*SuPIn
  1285. ENDIF
  1286. C*DBG-F call prfaible(ipt1,iel,nPIn0,mfaible)
  1287. *
  1288. 520 CONTINUE
  1289. *
  1290. * Ajustement au plus juste puis desactivation des segments lies
  1291. * a la rigidite du type 0
  1292. if (nelri0.ne.nelrig) then
  1293. nelrig = nelri0
  1294. segadj,xmatri
  1295. nbelem = nelri0
  1296. nbnn = 7
  1297. segadj,meleme
  1298. endif
  1299. segdes,xmatri
  1300. *
  1301. c write(ioimp,*) ' nb relation type 0 ',nelri0
  1302. * if (nelri0.eq.0) irigel(6,nrigel)=0
  1303. *
  1304. * Destruction des segments locaux devenus inutiles
  1305. segsup,mfaible
  1306. *
  1307. C* GOTO 1000
  1308.  
  1309. *=======================================================================
  1310. *= Fin du sous-programme :
  1311. *=======================================================================
  1312. 1000 CONTINUE
  1313. *
  1314. * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7
  1315. * puis desactivation du chpoint
  1316. if (vpocha(/1).ne.nelch) then
  1317. n = nelch
  1318. nc=vpocha(/2)
  1319. segadj,mpoval
  1320. nbnn = 1
  1321. nbelem = nelch
  1322. nbsous = 0
  1323. nbref = 0
  1324. segadj,ipt7
  1325. endif
  1326. * Desctivation de la matrice de raideur de contact
  1327. segdes,mrigid
  1328. *
  1329. * Reunion des relations portant sur le meme multiplicateur de lagrange
  1330. *
  1331. call impofu(MRIGID,MCHPOI)
  1332. *
  1333. * Nettoyage eventuel des termes petits dans les relations
  1334. *
  1335. * A voir plus tard. Frig3c suppose que la relation est complete.
  1336. *
  1337. ** call relasi(mrigid)
  1338. *
  1339. * Ecriture des resultats de l'operateur :
  1340. call ecrobj('RIGIDITE',mrigid)
  1341. call actobj('CHPOINT ',mchpoi,1)
  1342. call ecrobj('CHPOINT ',mchpoi)
  1343. *
  1344. 900 CONTINUE
  1345. end
  1346.  
  1347.  
  1348.  
  1349.  
  1350.  
  1351.  
  1352.  
  1353.  
  1354.  
  1355.  
  1356.  

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