Télécharger impo32.eso

Retour à la liste

Numérotation des lignes :

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

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